

!********************************************************************************************
!       THIS PROGRAM ACCOMPANIES Fried et al (2024)
!
!       "Understanding the Inequality and Welfare Impacts of Carbon Tax Policies"
!
!   This code is used to solve the transition from the steady state without the carbon tax to a new one with a tax
!   The user must solve for beginning and ending steady state prior to running program.
!   The use defines how the carbon tax revenue is used over the transition (e.g. how quickly the tax ramps up and how it is used over transition)
!   If you use this program, I would kindly ask that you cite my paper.
!*********************************************************************************************

PROGRAM hetero_ss
USE functions_mod
USE grid_interpolation
USE lambda0_optimizer_mod
USE retired_valuer_mod
USE working_valuer_mod
USE parameters_mod
USE parameters_energy_mod
USE NEQNF_INT

IMPLICIT NONE
INTERFACE
    function randu(iseed)
        INTEGER, intent(inout)   ::iseed
    DOUBLE PRECISION   ::randu
end function randu
END INTERFACE
include 'mpif.h'

parameter send_data_tag = 2001, return_data_tag = 2002

INTEGER            :: con_pos,kstart, capital_positive, tax_loop, non_clear(yr_T),temp,kit,i,s,it,it1,it4,it_tax
INTEGER            :: ac,out_loop,e_it,rebate_loop
INTEGER            :: sim_shocks(J,simulationsN,yr_T),sim_types(J,simulationsN,yr_T)
INTEGER            :: krow
INTEGER            :: shocks_start1,shocks_start2,routine
INTEGER            :: num_procs_run,seed,it2,it3,extra,optimal_baseline

INTEGER            :: my_id, root_process, ierr, status(MPI_STATUS_SIZE)
INTEGER            :: num_procs, an_id,vary_capital,update_lambda0_ig
INTEGER            :: num_procs_run2,share_it,share_it2,share_it3,yr
INTEGER            :: start_row,end_row,increment,large_inc,increment2,large_inc2
INTEGER            :: bad_combination_inner_nostop,bad_combination_outer,bad_combination_inner_stop


DOUBLE PRECISION   :: b,pricefactor(yr_T),xguess(3),solver_out(3),rate_use,ypp_temp
DOUBLE PRECISION   :: curv,diff_parameters(yr_T),diff_parameter,diff_parameter_v(6,yr_T),survive(J), dying(J),types(typesN)
DOUBLE PRECISION   :: growth,  pop_temp1(J),hours_jnr(yr_T),sim_earnings(J,simulationsN,yr_T)


DOUBLE PRECISION   :: v_end(typesN,shocksN,kgridN,J+1),v_retired_end(typesN,shocksN,kgridN,J+1)
DOUBLE PRECISION   :: agg_labor(J,typesN,yr_T),agg_capital(J+1,typesN,yr_T+1),agg_consumption(J,typesN,yr_T),agg_energy(J,typesN,yr_T)
DOUBLE PRECISION   :: kprime_dummy,labor_dummy

DOUBLE PRECISION   :: taxes_clear(max_outer_loops+1,yr_T),dummy,sim_v(simulationsN,yr_T),sim_total_income(J,simulationsN,yr_T),sim_after_tax_income(J,simulationsN,yr_T)


DOUBLE PRECISION   :: v_in(typesN,shocksN,kgridN),aggregates_start(31),aggregates_end(135),xguesses_check(3),check_dummy,check_dummy_out(3)
DOUBLE PRECISION   :: ctot80(kgridN,1),e80(kgridN,1),c80(kgridN,1), con(kgridN),con_dummy
DOUBLE PRECISION   :: sim_labor_tax_new(J,simulationsN,yr_T),sim_labor_tax_old(J,simulationsN,yr_T)
DOUBLE PRECISION   :: ss(typesN,(shocksN)/2,yr_T),ss1(typesN,(shocksN)/2,yr_T),ss_start(typesN,(shocksN)/2),ss_end(typesN,(shocksN)/2)

DOUBLE PRECISION   :: energy_miss,ind_eng_tax(yr_T),lump_sum_ig(yr_T)
DOUBLE PRECISION   :: c_upper(kgridN),avg_earnings1(yr_T),tax_dummy,labor_tax_dummy,capital_tax_dummy
DOUBLE PRECISION   :: increment_real,increment_real2

DOUBLE PRECISION   :: Ce_for_share_start,C_for_share_start,pe_start,&
r_start,lambdak_start

DOUBLE PRECISION   :: rebate_level_ig(yr_T),rebate_level1(yr_T),lambda1_benchmark

double precision aggregates(135,yr_T),parameters(yr_T,3)


DOUBLE PRECISION   :: ce_end(yr_T),C_end,Ce_start,C_start
DOUBLE PRECISION   :: lambda0(yr_T),lambda1(yr_T),lambda2(yr_T),lambdak(yr_T)

DOUBLE PRECISION   :: avg_earnings_ss_start(typesN,(shocksN)/2),avg_earnings_ss_end(typesN,(shocksN)/2)
DOUBLE PRECISION   :: K_start,K_end,Ec_start,Ec_end,tr_start,tr_end,N_start,N_end

DOUBLE PRECISION   :: hbar,partprod, K_ig(yr_T+1),E(yr_T),Ec(yr_T),Ec1(yr_T),Ep(yr_T), Ec_ig(yr_T),N_ig(yr_T), tr_ig(yr_T+1)
DOUBLE PRECISION   :: K(yr_T+1), N(yr_T), N1(yr_T),N2(yr_T),K1(yr_T),K2(yr_T),Ae,Atfp,sigmaa,gov_spend(yr_T)
DOUBLE PRECISION   :: K1_prime(yr_T+1),N1_prime(yr_T),Y(yr_T), w(yr_T), r(yr_T),kgridl, kgridu,  tr1(yr_T+1),sigmav,sigmaeps,rho
DOUBLE PRECISION   :: a,a1,tax_diff_one,tax_diff(yr_T),tax_diff_keep(yr_T),govt_taxes(yr_T),govt_taxes1(yr_T),govt_taxes2(yr_T),&
    govt_taxes3(yr_T),govt_taxes4(yr_T),tol,pop_start,tr(yr_T+1)
DOUBLE PRECISION   :: shocks_temp((shocksN)/2),stat_dist_temp((shocksN)/2),shock_pi_temp((shocksN)/2,(shocksN)/2)
DOUBLE PRECISION   :: shocks(shocksN),shock_pi(shocksN,shocksN,Jnr-1,typesN),shock_cumpi(shocksN,shocksN,J,typesN)
DOUBLE PRECISION   :: avg_earnings(yr_T),prices(11)
DOUBLE PRECISION   :: random(J,simulationsN)
DOUBLE PRECISION   :: avg_earnings_ss1(typesN,(shocksN)/2,yr_T),&
    avg_earnings_ss_ig(typesN,(shocksN)/2,yr_T),avg_earnings_ss(typesN,(shocksN)/2,yr_T)
DOUBLE PRECISION   :: ss_count(typesN,(shocksN)/2,yr_T)

DOUBLE PRECISION   :: sim_capital_start(J+1,simulationsN),sim_consumption_start(J,simulationsN),sim_labor_start(J,simulationsN),sim_energy_start(J,simulationsN),&
    sim_labor_earnings_start(J,simulationsN),sim_labor_tax_start(J,simulationsN),sim_capital_earnings_start(J,simulationsN),sim_capital_tax_start(J,simulationsN),&
    sim_total_after_tax_earnings_start(J,simulationsN),sim_lifetime_after_tax_income_start(simulationsN)



DOUBLE PRECISION, allocatable, dimension(:,:,:) :: sim_labor,sim_capital,sim_consumption,sim_energy,sim_rebate
DOUBLE PRECISION, allocatable, dimension(:,:,:,:,:):: v, v_retired
DOUBLE PRECISION, allocatable, dimension(:,:) :: zeros, population,dummy_zeros,sim_v_alloc
DOUBLE PRECISION, allocatable, dimension(:,:,:) :: sim_labor_alloc,sim_capital_alloc
DOUBLE PRECISION, allocatable, dimension(:,:,:) :: sim_consumption_alloc,sim_earnings_alloc,sim_energy_alloc,&
    sim_total_income_alloc,sim_after_tax_income_alloc
DOUBLE PRECISION, allocatable, dimension(:) :: zeros2,share_guesses,share_guesses2,share_guesses3
DOUBLE PRECISION, allocatable, dimension(:,:,:) :: v_alloc


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! MPI Parameters

root_process = 0
call MPI_INIT(ierr)
call MPI_COMM_RANK (MPI_COMM_WORLD, my_id, ierr)
call MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr)
num_procs_run=num_procs
num_procs_run2=num_procs
!Increments for simulations
if (simulationsN<(num_procs-1)) then
    num_procs_run2=simulationsN+1
end if
increment_real2=simulationsN/(num_procs_run2-1.0D0)
increment2=floor(increment_real2)
large_inc2=anint((increment_real2-increment2)*(num_procs_run2-1))

!Increments for value function calculation

if (kgridN<(num_procs-1)) then
    num_procs_run=kgridN+1
end if
increment_real=kgridN/(num_procs_run-1.0D0)
increment=floor(increment_real)
large_inc=anint((increment_real-increment)*(num_procs_run-1))


CALL ERSET(4,1,0)
CALL ERSET(3,1,0)
CALL ERSET(2,1,0)

household=0
allocate(v(typesN,shocksN,kgridN,Jnr,yr_T+1))
allocate(v_retired(typesN,shocksN,kgridN,J+1-jretired,yr_T+1))
allocate(sim_labor(J,simulationsN,yr_T))
allocate(sim_capital(J+1,simulationsN,yr_T+1))
allocate(sim_consumption(J,simulationsN,yr_T))
allocate(sim_energy(J,simulationsN,yr_T))
allocate(sim_rebate(J,simulationsN,yr_T))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Computational Parameters
a=0.07D0
a1=0.07D0
tol=0.00025D0
allocate(zeros2(max_outer_loops+1))
zeros2=0.0D0
do yr=1,yr_T
taxes_clear(:,yr)=zeros2
end do
deallocate(zeros2)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Utility Parameters
sigma1=2.0D0
sigma2=0.5D0

if(household==0) then
    !Individual
    gamma1=0.99074D0
!    ebar=0.001324D0*2.0D0
    ebar=0.001324D0
!    ebar=0.00D0
    beta=0.99502D0
    chi=73.3
    sigma2=0.5D0
    kgridl=-0.157D0
elseif (household==1) then
    !Household
    gamma1=0.99275D0
    ebar=0.0068D0
    beta=1.0085D0
    chi=90.0D0
    sigma2=0.55D0
    kgridl=0.0D0
else
    print*,"Set Household or Individual"
end if


!!!Set this to 1 if using unconstrained optimal as the baseline; 0 otherwise
optimal_baseline=0


rebate_to_old=1
negative_taxes=0
!Can rebate lead to negative effective taxes
rebate_negative_tax=1
!!!! 1 labor only, 2 capital only, 3 all
labor_capital_rebate_base=3


!No carbon tax but rebate revenue
no_carbon_tax=0


!Make rebate level grid capital grid (1=added as rebate level,2=added as rebate slope)
vary_capital=0



no_increase=1
!Doesn't allow labor tax to be larger than baseline
!-1 no longer constrains labor tax to be lower than previous
extra=0
!have the extra combinations on top, make sure decleration is right above



routine=3
!1 lump_sum_rebate
!2 throw_away
!3 cap_tax_rebate
!4 labor_rebate
!5 Rebate level

!Updated lambda0 when new parameters for IG
update_lambda0_ig=1
if(routine.NE.4) then
    update_lambda0_ig=0
end if

!Create tax that is 0 up to progressivity, then old tax after (routine should be 3)
flat_tax=0
if(flat_tax==1 .and. routine .NE. 3) then
print*,"UHOHOHOHOH flat tax on "
end if






!!!optimal
!do it2=1,yr_T
!    rebate_level(it2)=0.00D0
!    lambdak(it2)=.36D0
!    rebate_slope(it2)=0.0D0
!    lambda1(it2)=0.2325D0
!end do

!combination (routine 3)
do it2=1,yr_T
    lambdak(it2)=.36D0
    rebate_slope(it2)=0.0D0
    lambda1(it2)=0.2325D0
end do


!Overshoot progressivity(routine 3)
rebate_level=0.00
!lambda1=0.2325D0
!lambda1(1)=0.26D0
lambda1(1)=0.18D0


do it2=2,10
    lambda1(it2)=lambda1(it2-1)-(lambda1(1)-lambda1(11))/9.0D0
!    lambda1(it2)=0.2325D0
end do
do it2=10,yr_T
    rebate_level(it2)=0.0D0
    lambdak(it2)=.3125D0
end do

do it2=2,10
    lambdak(it2)=lambdak(it2-1)-(lambdak(1)-lambdak(yr_T))/9.0D0
end do



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Productivity Parameters
!Individual
hbar=0.00D0
partprod=1.0D0

if(household==0) then
epsilon1(:,1) = (/1.0D0, 1.05037955556458D0, 1.10168574952967D0, 1.15381030928342D0, 1.20663608757434D0, 1.26003734304819D0, 1.3138800980071D0, 1.36802257372016D0, 1.42231570281972D0, 1.47660371749329D0, 1.53072481133652D0, 1.584511871879D0, 1.63779327994293D0, 1.6903937711557D0, 1.74213535412359D0, 1.79283827899664D0, 1.84232204942572D0, 1.89040647024333D0, 1.93691272260084D0, 1.9816644577762D0, 2.02448890043734D0, 2.06521795181501D0, 2.10368928301267D0, 2.13974740856428D0, 2.1732447303485D0, 2.20404254208235D0, 2.23201198484835D0, 2.25703494445744D0, 2.27900488191109D0, 2.29782758879663D0, 2.31342186012357D0, 2.32572007787805D0, 2.33466869942757D0, 2.34022864584003D0, 2.34237558617614D0, 2.34110011486113D0, 2.33640782032644D0, 2.32831924422039D0, 2.31686973160448D0, 2.30210917366396D0, 2.28410164555347D0, 2.26292494305636D0, 2.23867002274592D0, 2.21144035128553D0, 2.18135117038086D0, 2.14852868468884D0/)
elseif(household==1) then
epsilon1(:,1) = (/1.0D0, 1.11972993442023D0, 1.23163188204661D0, 1.33628274802356D0, 1.43423225535758D0, 1.52600294491726D0, 1.61209017543323D0, 1.69296212349824D0, 1.76905978356708D0, 1.84079696795664D0, 1.90856030684589D0, 1.97270924827585D0, 2.03357605814965D0, 2.09146582023247D0, 2.14665643615159D0, 2.19939862539634D0, 2.24991592531814D0, 2.29840469113051D0, 2.34503409590901D0, 2.3899461305913D0, 2.4332556039771D0, 2.47505014272822D0, 2.51539019136855D0, 2.55430901228404D0, 2.59181268572273D0, 2.62788010979474D0, 2.66246300047226D0, 2.69548589158955D0, 2.72684613484296D0, 2.75641389979091D0, 2.7840321738539D0, 2.8095167623145D0, 2.83265628831737D0, 2.85321219286924D0, 2.87091873483891D0, 2.88548299095726D0, 2.89658485581726D0, 2.90387704187394D0, 2.90698507944442D0, 2.90550731670788D0, 2.89901491970559D0, 2.88705187234091D0, 2.86913497637924D0, 2.84475385144809D0, 2.81337093503703D0, 2.77442148249772D0/)
else
end if
!Firm
alpha1 = 0.3D0
alpha2 = 1.0D0-0.403D0
theta=0.03D0
delta = 0.079062D0
Atfp=1.0D0
Ae=Atfp*1.0D0

!taue=pe*0.326144D0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Heterogeneity Parameters
if(household==0) then
    rho=0.958D0
    sigmaa=0.065D0**0.5D0
    sigmav=(0.017D0/(1.0D0-rho**2.0D0))**0.5D0
    sigmaeps=0.081D0**0.5D0
elseif(household==1) then
    rho=0.972D0
    sigmaa=0.022D0**0.5D0
    sigmav=(0.012D0/(1.0D0-rho**2.0D0))**0.5D0
    sigmaeps=0.093D0**0.5D0
else
    print*,"Set Household"
end if



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Calculate Heterogeneity
do i =1,typesN
    types(i)=exp(-(sigmaa)+(real(i-1)/real(typesN-1))*2.0D0*sigmaa)
end do

! Stochastic parameters
CALL rouwenhorst(rho,sigmav,shock_pi_temp,shocks_temp,stat_dist_temp,(shocksN)/2)


shocks_temp=exp(shocks_temp)
shocks_temp=shocks_temp/(sum(stat_dist_temp*shocks_temp))

shocks(1:int(shocksN/2))=shocks_temp*exp(-sigmaeps)
shocks(int(shocksN/2)+1:int(shocksN))=shocks_temp*exp(sigmaeps)

if(my_id==1) then
print*,"shocks",shocks
end if
shocks_start1=int((((shocksN)/2)+1)/2)
shocks_start2=int(shocks_start1+((shocksN)/2))

!E to E
do it1=1,Jnr-1
do it2=1,typesN
shock_pi(1:int((shocksN)/2),1:int((shocksN)/2),it1,it2)=shock_pi_temp/2.0D0
shock_pi(int((shocksN)/2+1):(shocksN),1:int((shocksN)/2),it1,it2)=shock_pi_temp/2.0D0
shock_pi(1:int((shocksN)/2),int((shocksN)/2+1):(shocksN),it1,it2)=shock_pi_temp/2.0D0
shock_pi(int((shocksN)/2+1):(shocksN),int((shocksN)/2+1):(shocksN),it1,it2)=shock_pi_temp/2.0D0
end do
end do





!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Population Parameters
growth = 0.011D0
if(household==0) then
pop_temp1 =(/197316.D0, 197141.D0, 196959.D0, 196770.D0, 196580.D0, 196392.D0, 196205.D0, 196019.D0, 195830.D0, 195634.D0, 195429.D0, 195211.D0, 194982.D0, 194739.D0, 194482.D0, 194211.D0, 193924.D0, 193619.D0, 193294.D0, 192945.D0, 192571.D0, 192169.D0, 191736.D0, 191271.D0, 190774.D0, 190243.D0, 189673.D0, 189060.D0, 188402.D0, 187699.D0, 186944.D0, 186133.D0, 185258.D0, 184313.D0, 183290.D0, 182181.D0, 180976.D0, 179665.D0, 178238.D0, 176689.D0, 175009.D0, 173187.D0, 171214.D0, 169064.D0, 166714.D0, 164147.D0, 161343.D0, 158304.D0, 155048.D0, 151604.D0, 147990.D0, 144189.D0, 140180.D0, 135960.D0, 131532.D0, 126888.D0, 122012.D0, 116888.D0, 111506.D0, 105861.D0, 99957.D0, 93806.D0, 87434.D0, 80882.D0, 74204.D0, 67462.D0, 60721.D0, 54053.D0,	47533.D0,	41241.D0, 35259.D0, 29663.D0, 24522.D0, 19890.D0, 15805.D0, 12284.D0, 9331.D0, 6924.D0, 5016.D0, 3550.D0, 2454.D0 /)
do i=1,(J-1)
   survive(i) = pop_temp1(i+1)/pop_temp1(i)
end do
survive(J) = 0.0D0
adult_equivalence=1.0D0
else
survive=(/1.0D0, 0.999443429520098D0, 0.999453048952094D0, 0.999474759951225D0, 0.999507517658172D0, 0.999550180294631D0, 0.99958688938571D0, 0.999616966794992D0, 0.999637259341625D0, 0.999646435563957D0, 0.999642137784891D0, 0.999638705313936D0, 0.999628264835928D0, 0.999611266702816D0, 0.999587623092204D0, 0.99955347402391D0, 0.999517484919676D0, 0.999474904683543D0, 0.999424216406889D0, 0.999368982231295D0, 0.999311637585569D0, 0.999247262201619D0, 0.99917475362204D0, 0.999096164795857D0, 0.999015366949988D0, 0.99892772761533D0, 0.998832359878012D0, 0.998729478390282D0, 0.998631176268071D0, 0.998535878125485D0, 0.998443264975948D0, 0.998338654631507D0, 0.998221906285934D0, 0.998078357864358D0, 0.997913954710353D0, 0.997725564043014D0, 0.997510605728885D0, 0.997277513805754D0, 0.99702940364917D0, 0.996758720419625D0, 0.996463887451456D0, 0.996136496024502D0, 0.995758248748826D0, 0.995343091935131D0, 0.994874730958277D0, 0.994354016314192D0, 0.993755337681475D0, 0.993065683870379D0, 0.992293047744931D0, 0.991433174915763D0, 0.990457448940283D0, 0.989303075213917D0, 0.987964810901374D0, 0.986475935412401D0, 0.984833077329247D0, 0.982973697702629D0, 0.980769120672249D0, 0.978166137633033D0, 0.975176536129091D0, 0.971754588528722D0, 0.967817497628447D0, 0.963253722024177D0, 0.957955278732916D0, 0.951933505633426D0, 0.945105096525417D0, 0.937431585348977D0, 0.92887640431547D0, 0.919444350437676D0, 0.909124748145856D0, 0.89794551325945D0, 0.885974030114678D0, 0.873304919941224D0, 0.859956433479889D0, 0.842665940409714D0, 0.824214308977887D0, 0.804792055183895D0, 0.785116734289799D0, 0.765756410717599D0, 0.747202819257837D0, 0.729458167443388D0, 0.713669733378468D0/)
survive(J) = 0.0D0
adult_equivalence=(/1.21998398833892D0, 1.26899136580154D0, 1.32056246723579D0, 1.3734138700101D0, 1.42641978111077D0, 1.47860215039009D0, 1.52912105574409D0, 1.57726536021961D0, 1.62244364105098D0, 1.66417539062616D0, 1.70208248938226D0, 1.73588095063074D0, 1.76537293731188D0, 1.79043905067891D0, 1.81103089091152D0, 1.82716388965886D0, 1.83891041451209D0, 1.84639314540638D0, 1.84977872295234D0, 1.84927166869702D0, 1.84510857731432D0, 1.83755258072496D0, 1.8268880841459D0, 1.81341577406915D0, 1.79744789817028D0, 1.77930381714618D0, 1.75930582848249D0, 1.73777526215036D0, 1.71502884823281D0, 1.69137535648057D0, 1.66711250779728D0, 1.64252415765435D0, 1.61787775143512D0, 1.5934220517087D0, 1.5693851374331D0, 1.54597267508807D0, 1.52336646173711D0, 1.50172324001931D0, 1.48117378507042D0, 1.46182226337349D0, 1.44374586353911D0, 1.42699469901482D0, 1.41159198272448D0, 1.39753447363654D0, 1.3847931952624D0, 1.37331442608375D0, 1.36302096190969D0, 1.35381365016324D0, 1.34557319609731D0, 1.33816224094031D0, 1.33142771197095D0, 1.32520344452275D0, 1.31931307591802D0, 1.31357321133145D0, 1.30779686158245D0, 1.30179715285813D0, 1.29539130836514D0, 1.28840490191053D0, 1.28067638341311D0, 1.27206187634362D0, 1.26244024709437D0, 1.25171844627878D0, 1.23983712196005D0, 1.22677650480938D0, 1.21256256519385D0, 1.19727344219358D0, 1.18104614454827D0, 1.16408352353373D0, 1.14666151776743D0, 1.12913666994348D0, 1.11195391549768D0, 1.09896503226485D0, 1.08988645284858D0, 1.0816144615019D0, 1.07430423674565D0, 1.0681188559846D0, 1.06322914430733D0, 1.0598135018135D0, 1.0580577094685D0, 1.05815471348597D0, 1.06030438823754D0/)
!adult_equivalence=1.0D0
end if

dying = 1.D0 - survive
allocate(population(J,2))
population=0.0D0
population(1,1) = 1.0D0
do i=2,(J)
   population(i,1) = survive(i-1)*population((i-1),1)/(1+growth)
end do
pop_start=population(1,1)/sum(population(1:J,1))
population(:,2) = population(:,1)/sum(population(:,1))



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Grids
kgrid(1)=kgridl
kgrid(kgridnegN+1)=0.0D0
curv = 1.5D0
do ac=2,kgridnegN+1
	kgrid(kgridnegN-ac+2)=kgrid(kgridnegN+1)+(kgridl)*((ac-1.0D0)/(kgridnegN))**curv
end do

kgridu = 25.0D0
curv = 1.35D0
do ac=kgridnegN+2,kgridN
	kgrid(ac)=kgrid(kgridnegN+1)+(kgridu-kgrid(kgridnegN+1))*((ac-kgridnegN-1.0D0)/(kgridN-kgridnegN-1.0D0))**curv
end do
kgrid(kgridN)=kgridu
capital_positive = minloc(kgrid,1, mask = kgrid.ge.0)
kstart=capital_positive



!kgridl = 0.0D0
!kgridu = 101.0D0
!kstart = 1
!!Curved grids to handle outliers
!curv = 2.25D0
!kgrid(1)=kgridl
!
!do ac=2,kgridN
!	kgrid(ac)=kgrid(1)+(kgridu-kgridl)*((ac-1.0D0)/(kgridN-1.0D0))**curv
!end do
!kgrid(kgridN)=kgridu


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Social Security Payment Information
b=0.45D0

!!!!!!!!!!!!!!!!!!!!!!!!!!!! Government Parameters
!gs=.157D0
govt_taxes=0.0D0
lambda1_benchmark=.031D0
!lambda1=.031D0

if(my_id==1) then
print*,"here"
end if
 OPEN(227, FILE="/<filepath starting steady state>/aggregates_output.dat",FORM="formatted")
    read(227,887) (aggregates_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steady state>/ss_output.dat",FORM="formatted")
    read(227,887) (ss_start)
    close(227)

 OPEN(227, FILE="/<filepath starting steady state>/avg_earnings_ss_output.dat",FORM="formatted")
    read(227,887) (avg_earnings_ss_start)
    close(227)

gov_spend=aggregates_start(7)
lambda0=aggregates_start(6)

Ce_for_share_start=aggregates_start(14)
C_for_share_start=aggregates_start(15)
pe_start=aggregates_start(16)
!    ss_start=aggregates_start(9)
!    ss_start(2)=aggregates_start(10)
avg_earnings=aggregates_start(17)
!initial guess
!ss_cap=avg_earnings*2.42D0
avg_earnings_scale=(aggregates_start(17))*(1.0D0-aggregates_start(8)/2.0D0)


K_start=aggregates_start(3)
N_start=aggregates_start(4)
Ec_start=aggregates_start(22)
tr_start=aggregates_start(5)



!!!!!!!!!!!!!!Read in ending
 OPEN(227, FILE="/<filepath ending steady state>/transition_steady_state/optimal/optimal.dat",FORM="formatted")
    read(227,887) (aggregates_end)
    close(227)

 OPEN(227, FILE="/<filepath ending steady state>/transition_steady_state/optimal/ss_output.dat",FORM="formatted")
    read(227,887) (ss_end)
    close(227)

 OPEN(227, FILE="/<filepath ending steady state>/transition_steady_state/optimal/avg_earnings_ss_output.dat",FORM="formatted")
    read(227,887) (avg_earnings_ss_end)
    close(227)

 OPEN(227, FILE="/<filepath ending steady state>/transition_steady_state/optimal/v_output.dat",FORM="formatted")
    read(227,887) (v_end)
    close(227)

 OPEN(227, FILE="/<filepath ending steady state>/transition_steady_state/optimal/v_retired_output.dat",FORM="formatted")
    read(227,887) (v_retired_end)
    close(227)

v_retired(:,:,:,:,yr_T+1)=v_retired_end(:,:,:,Jnr:J+1)
v(:,:,:,:,yr_T+1)=v_end(:,:,:,1:Jnr)
K_end=aggregates_end(3)
N_end=aggregates_end(4)
Ec_end=aggregates_end(79)
tr_end=aggregates_end(5)
if (routine==3) then
    lambda0=aggregates_end(13)

end if
!!!No tax
!taue=0.0D0

!!Standard tax
taue=pe_start*0.482D0

!Small tax
!taue=pe_start*0.482D0*0.75D0

!!large tax
!taue=pe_start*0.482D0*1.25D0
!!!! optimal labor: .792678 .121 capital: .116
!!!! benchmark labor: .826689 .031 capital: .36
lambda0_start=aggregates_start(6)
lambda1_start=aggregates_start(13)

rebate_level_ig=0.0D0
r_start=aggregates_start(2)
lambdak_start=aggregates_start(12)

if(no_carbon_tax==1) then
!taue=pe_start*0.482D0

print*,"NOT TAXING CARBON"
!!Base on st st
!gov_spend=aggregates_start(7)-(aggregates_start(27)+aggregates_start(22))*taue
!Base on optimal
if(optimal_baseline==0) then
gov_spend=aggregates_start(7)-(0.006164949D0+0.025624986D0)*taue
else
!    print*,"STOOOOOOOOOOOOOOPPPPP AND FILL IN FROM OPTIMAL"
gov_spend=aggregates_start(7)-(0.006167261D0+0.026070431D0)*taue
end if
taue=0.0D0
!gov_spend=aggregates_start(7)*(1.0D0-.0785)
end if

if(my_id==1) then
print*,"here2"
end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Initial Guesses



avg_earnings_ss_ig(:,:,1)=avg_earnings_ss_start
K_ig(1)=K_start
N_ig(1)=N_start
Ec_ig(1)=Ec_start
tr_ig(1)=tr_start

avg_earnings_ss_ig(:,:,yr_T)=avg_earnings_ss_end
K_ig(yr_T)=K_end
K_ig(yr_T+1)=K_end
N_ig(yr_T)=N_end
Ec_ig(yr_T)=Ec_end
tr_ig(yr_T)=tr_end
tr_ig(yr_T+1)=tr_end

do it=2,40
    avg_earnings_ss_ig(:,:,it)=avg_earnings_ss_ig(:,:,it-1)+(avg_earnings_ss_end-avg_earnings_ss_start)/(40.0D0-1.0D0)
    K_ig(it)=K_ig(it-1)+(K_end-K_start)/(40.0D0-1.0D0)
    N_ig(it)=N_ig(it-1)+(N_end-N_start)/(40.0D0-1.0D0)
    Ec_ig(it)=Ec_ig(it-1)+(Ec_end-Ec_start)/(40.0D0-1.0D0)
    tr_ig(it)=tr_ig(it-1)+(tr_end-tr_start)/(40.0D0-1.0D0)
end do

do it=41,yr_T
    avg_earnings_ss_ig(:,:,it)=avg_earnings_ss_ig(:,:,it-1)
    K_ig(it)=K_ig(it-1)
    N_ig(it)=N_ig(it-1)
    Ec_ig(it)=Ec_ig(it-1)
    tr_ig(it)=tr_ig(it-1)
end do


!Combo 20
K_ig=(/1.74962985238568D0, 1.74635380588454D0, 1.74387280034879D0, 1.74200337476381D0, 1.74064459067011D0, 1.73973392768733D0, 1.73917791020142D0, 1.73890175435927D0, 1.73882258096399D0, 1.73885812438711D0, 1.73892011906079D0, 1.7389962674026D0, 1.73907849620672D0, 1.73916147560151D0, 1.73924333651032D0, 1.73932225340548D0, 1.739400637724D0, 1.73947507843045D0, 1.73954573379387D0, 1.73961314790811D0, 1.73967843244128D0, 1.73974050872237D0, 1.73979863609769D0, 1.73985270859278D0, 1.73990402589078D0, 1.73995147430333D0, 1.73999575901219D0, 1.74003684212446D0, 1.74007522923563D0, 1.74010964043269D0, 1.74014094253004D0, 1.7401693219809D0, 1.74019483894151D0, 1.74021640140561D0, 1.74023528128989D0, 1.74025223352396D0, 1.74026674611273D0, 1.74027929372982D0, 1.74029016747183D0, 1.74029927333535D0, 1.74030782531089D0, 1.74031561130778D0, 1.74032257326055D0, 1.74032869844641D0, 1.74033424371008D0, 1.74033929951953D0, 1.74034400449626D0, 1.74034804117459D0, 1.74035126621562D0, 1.74035405670047D0, 1.74035630134697D0, 1.74035818228591D0, 1.74035969018454D0, 1.7403607268434D0, 1.74036139980618D0, 1.74036164151058D0, 1.74036162071287D0, 1.74036138518077D0, 1.74036111490223D0, 1.74036083191278D0, 1.74036058976154D0, 1.74036025738792D0, 1.74035999735903D0, 1.74035986197735D0, 1.74035980502428D0, 1.74035990628197D0, 1.74036009764462D0, 1.74036039359615D0, 1.74036078891368D0, 1.74036130324901D0, 1.74036189462346D0, 1.74036257420343D0, 1.74036330853666D0, 1.74036411785473D0, 1.74036496934706D0, 1.74036589084548D0, 1.74037308257444D0, 1.74037973619402D0, 1.74038729083456D0, 1.74039583946121D0, 1.74040557990017D0, 1.74041671053527D0, 1.74042929277028D0, 1.74044388037696D0, 1.74046009123355D0, 1.74027113991225D0, 1.74027551839842D0, 1.74027964087857D0, 1.74028348005984D0, 1.74028714683878D0, 1.74029068201804D0, 1.74029396929391D0, 1.74029709658113D0, 1.74029991791949D0, 1.74030241621651D0, 1.74030457360047D0, 1.74030638382679D0, 1.74030779865132D0, 1.74030879999506D0, 1.74030934965865D0, 1.7403094303477D0/)
N_ig=(/0.522071459258623D0, 0.522470600714396D0, 0.522796168069277D0, 0.523070667764561D0, 0.523309193077685D0, 0.523494195466028D0, 0.523641261826954D0, 0.523741747528235D0, 0.523798039774676D0, 0.523802934605059D0, 0.523803027445904D0, 0.523799203786302D0, 0.52379424282848D0, 0.523788793445268D0, 0.523783469572706D0, 0.523779838807748D0, 0.523774720524889D0, 0.523770626784097D0, 0.523767473204D0, 0.523766583995245D0, 0.523765452696817D0, 0.523763334788573D0, 0.523761114468545D0, 0.523760535744376D0, 0.523759440097926D0, 0.523759100065171D0, 0.523758601575757D0, 0.523758745831676D0, 0.523758072378523D0, 0.523757841896825D0, 0.523757955363598D0, 0.523758262836883D0, 0.523757790293987D0, 0.523758134787055D0, 0.523759114954153D0, 0.52375949892608D0, 0.523760435550756D0, 0.523761570916027D0, 0.523762753009229D0, 0.523763944059589D0, 0.523765167964602D0, 0.523766108931646D0, 0.523766898828546D0, 0.523767600416291D0, 0.523768275557926D0, 0.523768666840256D0, 0.523769009936974D0, 0.523769192948686D0, 0.523769220675712D0, 0.523769134592998D0, 0.523769032881408D0, 0.523768843837994D0, 0.523768637483839D0, 0.523768391901068D0, 0.523768115717229D0, 0.523767827809224D0, 0.523767569229292D0, 0.523767317855925D0, 0.523767101344339D0, 0.52376690714943D0, 0.523766741468271D0, 0.523766588787031D0, 0.523766465586116D0, 0.523766355337958D0, 0.523766260354775D0, 0.523766188817033D0, 0.523766126611475D0, 0.523766055535494D0, 0.523766029799234D0, 0.523766019788963D0, 0.523766021346853D0, 0.523766034359379D0, 0.523766058755709D0, 0.523766094058061D0, 0.523766139495318D0, 0.523766195588064D0, 0.523770300221745D0, 0.523770890394752D0, 0.523771525164274D0, 0.523772319022447D0, 0.523773255285367D0, 0.523774116300861D0, 0.523775513390772D0, 0.523776379542621D0, 0.523777895127374D0, 0.523772086096511D0, 0.523771956231104D0, 0.523771768196734D0, 0.523771762893691D0, 0.523771886835635D0, 0.523771812431899D0, 0.523771891520314D0, 0.523771762513386D0, 0.5237716095691D0, 0.523771440560533D0, 0.523771286949585D0, 0.523771064846824D0, 0.523770848118657D0, 0.523770586263069D0, 0.523770339100135D0/)
K_ig(1)=K_start
tr_ig(1)=tr_start
K_ig(yr_T)=K_end
K_ig(yr_T+1)=K_end
tr_ig(yr_T)=tr_end
tr_ig(yr_T+1)=tr_end
tr_ig=(/0.0231359281706988D0, 0.0230917072561257D0, 0.0230553595100321D0, 0.0230266246124503D0, 0.0230044872825005D0, 0.0229880332508093D0, 0.0229762632650157D0, 0.0229683083152848D0, 0.0229632776516319D0, 0.0229603452357237D0, 0.0229587002733272D0, 0.0229573146967692D0, 0.0229561573129817D0, 0.0229551400727629D0, 0.0229542436139538D0, 0.0229534251875956D0, 0.0229527571936176D0, 0.022952198456391D0, 0.0229517434881992D0, 0.0229514168662255D0, 0.0229512672606929D0, 0.0229512770254652D0, 0.0229514702808213D0, 0.0229518695460884D0, 0.0229524866371736D0, 0.0229532996033243D0, 0.0229543058158036D0, 0.0229554878660379D0, 0.0229568180735791D0, 0.0229582431631624D0, 0.0229597500573954D0, 0.0229612865634374D0, 0.0229628270115006D0, 0.0229643434286219D0, 0.0229658303223612D0, 0.0229672888392024D0, 0.0229686825227386D0, 0.0229700033227413D0, 0.0229712446978757D0, 0.0229723932374996D0, 0.0229734810589209D0, 0.0229744962390735D0, 0.0229754393829262D0, 0.0229762987758586D0, 0.022977089410216D0, 0.0229778011361921D0, 0.0229784539275852D0, 0.0229790388036239D0, 0.0229795464727822D0, 0.0229800030765448D0, 0.0229804017045064D0, 0.0229807483947508D0, 0.0229810415775698D0, 0.0229812685487542D0, 0.0229814406460676D0, 0.0229815521422788D0, 0.0229816091106732D0, 0.0229816154545217D0, 0.0229815786409127D0, 0.0229815098461512D0, 0.0229814197737469D0, 0.0229813033006509D0, 0.0229811746151754D0, 0.0229810452930031D0, 0.0229809174099991D0, 0.0229808015250246D0, 0.0229806926699926D0, 0.0229805974033545D0, 0.022980517588883D0, 0.0229804555474072D0, 0.0229804055577395D0, 0.0229803705493156D0, 0.0229803460545792D0, 0.0229803346444214D0, 0.0229803303118034D0, 0.0229803342146078D0, 0.0229799127912039D0, 0.02297993870126D0, 0.0229799725766279D0, 0.0229800151017732D0, 0.0229800672136663D0, 0.0229801299473987D0, 0.022980203329273D0, 0.0229802903320685D0, 0.0229776285334239D0, 0.022977681840942D0, 0.0229777755878772D0, 0.0229778712032808D0, 0.022977968638393D0, 0.0229780685058587D0, 0.0229781707127269D0, 0.0229782753234545D0, 0.0229783827211495D0, 0.0229784925592992D0, 0.0229786045509396D0, 0.0229787186551204D0, 0.0229788350151906D0, 0.0229789531329176D0, 0.0229790727279668D0, 0.0229791932733051D0, 0.0229793135696444D0/)
Ec_ig=(/0.0062001983392168D0, 0.0061919916150009D0, 0.0061850048064834D0, 0.00617940655382D0, 0.0061748498457202D0, 0.0061713072641087D0, 0.0061686807393393D0, 0.0061668575391072D0, 0.0061658346382574D0, 0.0061655111172587D0, 0.0061653018549373D0, 0.0061651974404289D0, 0.0061651438252957D0, 0.00616511528734D0, 0.0061651037677593D0, 0.0061650761828229D0, 0.006165074844917D0, 0.0061650820856669D0, 0.0061650814601024D0, 0.0061650844046266D0, 0.0061650973211606D0, 0.0061651072660241D0, 0.0061651139266215D0, 0.0061651180671816D0, 0.0061651235889656D0, 0.0061651236875835D0, 0.0061651186022318D0, 0.0061651154552854D0, 0.0061651164450807D0, 0.0061651077971808D0, 0.0061650979758474D0, 0.0061650853249825D0, 0.0061650745542045D0, 0.0061650544313825D0, 0.0061650261665036D0, 0.0061649983572751D0, 0.0061649676631256D0, 0.0061649329136483D0, 0.0061648981986119D0, 0.0061648498346399D0, 0.0061648764104341D0, 0.0061649015259025D0, 0.0061649236455905D0, 0.0061649411103007D0, 0.0061649569549944D0, 0.0061649687745057D0, 0.0061649842997759D0, 0.0061649988922164D0, 0.0061650069633305D0, 0.0061650153447666D0, 0.0061650200891135D0, 0.0061650245517696D0, 0.0061650295817403D0, 0.006165032672427D0, 0.0061650360303449D0, 0.0061650363895192D0, 0.0061650370150633D0, 0.006165034984766D0, 0.0061650328832663D0, 0.0061650305529354D0, 0.006165029947078D0, 0.0061650276726787D0, 0.0061650250926935D0, 0.0061650229527252D0, 0.0061650203930099D0, 0.0061650189045985D0, 0.0061650172238195D0, 0.006165015925641D0, 0.0061650147925586D0, 0.006165014324254D0, 0.0061650138404149D0, 0.0061650139766906D0, 0.0061650138971255D0, 0.0061650145341398D0, 0.0061650149020915D0, 0.0061650154095238D0, 0.006165035277804D0, 0.0061650338768041D0, 0.0061650324326094D0, 0.0061650306082837D0, 0.0061650284041236D0, 0.0061650261563172D0, 0.006165022250582D0, 0.0061650176839224D0, 0.0061650109983659D0, 0.006164959471215D0, 0.0061649623305893D0, 0.0061649651532346D0, 0.0061649676617247D0, 0.0061649705555735D0, 0.0061649731782749D0, 0.0061649758607899D0, 0.0061649786825569D0, 0.0061649814697259D0, 0.0061649841159146D0, 0.0061649868896254D0, 0.0061649895244511D0, 0.0061649921127339D0, 0.006164994765072D0, 0.0061649974973842D0/)
avg_earnings=(/0.554670990626827D0, 0.55463704469159D0, 0.554627675214206D0, 0.554641805270613D0, 0.554680449401954D0, 0.554724560856888D0, 0.554776442378509D0, 0.554822406792864D0, 0.554855611271283D0, 0.554862846079561D0, 0.554869153142998D0, 0.554873976819386D0, 0.554878590182583D0, 0.554882930275334D0, 0.55488742447029D0, 0.554892804876162D0, 0.554896862000307D0, 0.554901283677141D0, 0.554905990062393D0, 0.55491230295974D0, 0.554918155616964D0, 0.554922821924147D0, 0.554927047745718D0, 0.554932088333049D0, 0.554936434139957D0, 0.55494102622427D0, 0.554945087845636D0, 0.55494935793519D0, 0.55495272909849D0, 0.554956037686271D0, 0.554959289990137D0, 0.554962386785516D0, 0.554964612846019D0, 0.554967032571334D0, 0.554969751811257D0, 0.554971706570403D0, 0.554973829960956D0, 0.554975893984088D0, 0.554977811463876D0, 0.554979551735626D0, 0.554981087845528D0, 0.554982364056723D0, 0.554983469306951D0, 0.55498445258523D0, 0.554985379814918D0, 0.55498606302339D0, 0.554986673004402D0, 0.554987126996169D0, 0.554987402490171D0, 0.55498755787701D0, 0.554987665147838D0, 0.554987660140525D0, 0.554987623550127D0, 0.554987508774737D0, 0.554987348896273D0, 0.554987141202528D0, 0.554986932612125D0, 0.554986709609404D0, 0.554986512141189D0, 0.554986331101813D0, 0.554986176927178D0, 0.554986023127878D0, 0.554985896346472D0, 0.554985794997778D0, 0.554985717677483D0, 0.554985675581304D0, 0.554985649943911D0, 0.554985622878173D0, 0.554985644793905D0, 0.554985690247428D0, 0.554985752105275D0, 0.55498583120369D0, 0.554985924463668D0, 0.55498603326285D0, 0.554986153628033D0, 0.554986288642702D0, 0.554991177028231D0, 0.554992162742403D0, 0.554993270858227D0, 0.554994585825725D0, 0.554996110681474D0, 0.554997739683027D0, 0.55499989308783D0, 0.555001938627464D0, 0.555004460005632D0, 0.554988219015514D0, 0.554988206046215D0, 0.554988125056873D0, 0.554988220322763D0, 0.554988439768544D0, 0.554988454509795D0, 0.554988602981712D0, 0.554988551527879D0, 0.554988471857127D0, 0.554988366389569D0, 0.554988262670449D0, 0.554988085610948D0, 0.55498790302691D0, 0.554987666468969D0, 0.554987426828813D0/)

do it=50,yr_T
    avg_earnings_ss_ig(:,:,it)=avg_earnings_ss_end
end do

avg_earnings_ss_ig(1,1,:)=(/0.318969145548654D0, 0.318939578746734D0, 0.318973947493249D0, 0.319043088464792D0, 0.319112628480207D0, 0.319197089610624D0, 0.319283190510576D0, 0.319362562240597D0, 0.31943343505495D0, 0.31947234088151D0, 0.319511488020737D0, 0.319547985432577D0, 0.319572254711443D0, 0.319598708644628D0, 0.319621845127409D0, 0.319647201398634D0, 0.31966378632131D0, 0.319682317904486D0, 0.319706117190733D0, 0.319714648739128D0, 0.319732839274349D0, 0.319751795746664D0, 0.319757834534176D0, 0.319766739968258D0, 0.319777396390117D0, 0.319791789264448D0, 0.319802676137458D0, 0.319814194462568D0, 0.319820824297454D0, 0.31982618830421D0, 0.319835464861491D0, 0.31983930170957D0, 0.319843137343838D0, 0.319845266379299D0, 0.319848447509991D0, 0.319851671310167D0, 0.319851875754919D0, 0.319854826234742D0, 0.319857287794262D0, 0.319856210940106D0, 0.319857425606724D0, 0.31985814199121D0, 0.319858687687791D0, 0.319860532636502D0, 0.319861386652599D0, 0.31986148314693D0, 0.319861808378199D0, 0.319861917339233D0, 0.319861810274899D0, 0.319861744100383D0, 0.319861725254662D0, 0.319861596698968D0, 0.319861567999363D0, 0.319861422700178D0, 0.319861268670951D0, 0.319861116487022D0, 0.319860965437D0, 0.319860822194818D0, 0.319860692798031D0, 0.319860579516654D0, 0.319860482797583D0, 0.319860391472322D0, 0.319860311165097D0, 0.319860252401297D0, 0.319860212497643D0, 0.319860192134036D0, 0.319860186609657D0, 0.319860195934671D0, 0.319860213486528D0, 0.319860243385625D0, 0.319860283155221D0, 0.319860333555537D0, 0.319860391490273D0, 0.319860457909887D0, 0.319860530614607D0, 0.31986061278711D0, 0.319860704508606D0, 0.319860805417123D0, 0.319860916480432D0, 0.319861038572312D0, 0.319861172612385D0, 0.319861322712069D0, 0.319861489229083D0, 0.319861675958991D0, 0.319861881905919D0, 0.319862113462796D0, 0.319862375514586D0, 0.319862663144059D0, 0.319862997113719D0, 0.319863373458288D0, 0.319863799559587D0, 0.319864279060665D0, 0.319864821511561D0, 0.319865451149152D0, 0.31986614895791D0, 0.319866950682453D0, 0.319867870772199D0, 0.319868808366263D0, 0.319869999230175D0, 0.319871408114308D0/)
avg_earnings_ss_ig(2,1,:)=(/0.486469170003321D0, 0.48648727072704D0, 0.486506582666342D0, 0.486531432692842D0, 0.486593819764119D0, 0.486650207112261D0, 0.486717539206218D0, 0.486781827303022D0, 0.486832882618461D0, 0.48685892672965D0, 0.486865503194539D0, 0.486879104329951D0, 0.486888291366016D0, 0.486903117136368D0, 0.486905314946682D0, 0.486909069248578D0, 0.486912613910106D0, 0.486917622635516D0, 0.486923431376465D0, 0.486932898625364D0, 0.486945641797169D0, 0.486948015333918D0, 0.48695069250115D0, 0.486953730474159D0, 0.486952137951079D0, 0.48695576832876D0, 0.486958455627434D0, 0.486960288108922D0, 0.486960274615477D0, 0.486959852370454D0, 0.486960324152229D0, 0.486962371049925D0, 0.48696213693453D0, 0.486963060314378D0, 0.486963850465923D0, 0.48696611166119D0, 0.486965341014534D0, 0.486968528453912D0, 0.486968215222686D0, 0.486966708951146D0, 0.486966897477599D0, 0.486966966367578D0, 0.486969610152629D0, 0.486970233292273D0, 0.486971248449956D0, 0.486971705013936D0, 0.48697134585124D0, 0.486971262008956D0, 0.486971312121754D0, 0.48697126479622D0, 0.486971091205179D0, 0.486970488456903D0, 0.486970322657306D0, 0.486970155003017D0, 0.486969979258456D0, 0.486969794325391D0, 0.486969596332655D0, 0.486969407225168D0, 0.486969233902106D0, 0.486969077175679D0, 0.486968943499088D0, 0.486968817290728D0, 0.486968701697349D0, 0.48696861800603D0, 0.486968559193885D0, 0.486968529386913D0, 0.486968520751262D0, 0.486968532918234D0, 0.486968559037333D0, 0.486968604579044D0, 0.486968665124053D0, 0.486968740353308D0, 0.486968826810182D0, 0.48696892735093D0, 0.486969038313106D0, 0.48696916171164D0, 0.486969299459511D0, 0.48696944954887D0, 0.486969614749765D0, 0.486969793875019D0, 0.48696999360869D0, 0.486970214965057D0, 0.486970461601497D0, 0.4869707409462D0, 0.486971049753247D0, 0.486971395318514D0, 0.486971784062325D0, 0.486972212903304D0, 0.486972706736479D0, 0.486973259149818D0, 0.486973880597092D0, 0.48697458607209D0, 0.48697539266337D0, 0.486976311476624D0, 0.486977337055935D0, 0.486978525810808D0, 0.486979906478807D0, 0.486981355592821D0, 0.486983237156368D0, 0.486985413452259D0/)
avg_earnings_ss_ig(1,2,:)=(/0.371327141111113D0, 0.371305726109451D0, 0.37133473235087D0, 0.371391900564438D0, 0.371472535380909D0, 0.371552826062144D0, 0.37163106346846D0, 0.371705788799077D0, 0.371773247184842D0, 0.371819731051895D0, 0.371860347059456D0, 0.371889106193745D0, 0.37191791449244D0, 0.371939470537037D0, 0.371962664256641D0, 0.371982671489642D0, 0.372001090390944D0, 0.372019285428441D0, 0.372037700858975D0, 0.372058733978379D0, 0.372073742654701D0, 0.372087079642399D0, 0.372100719293716D0, 0.372113592487431D0, 0.372124752515525D0, 0.372134276275282D0, 0.372143185472392D0, 0.372151037677195D0, 0.372157519718839D0, 0.372164710545608D0, 0.37217168620235D0, 0.372177530924339D0, 0.372179939198312D0, 0.372183952015349D0, 0.372187966017113D0, 0.372190624655619D0, 0.37219287434673D0, 0.372195151469308D0, 0.372197473595829D0, 0.372199892849966D0, 0.372201093405061D0, 0.372202775362903D0, 0.372203424673111D0, 0.372205240809284D0, 0.372206761018612D0, 0.372207275695811D0, 0.372208054211841D0, 0.372208525538838D0, 0.372208846870918D0, 0.37220894495332D0, 0.372208964282817D0, 0.372209009732568D0, 0.372208964718633D0, 0.37220883878186D0, 0.372208687453054D0, 0.37220851052973D0, 0.372208333304096D0, 0.372208161528473D0, 0.372208004705343D0, 0.372207864547537D0, 0.372207743487666D0, 0.372207624959968D0, 0.372207517373715D0, 0.372207434027479D0, 0.372207363892434D0, 0.372207315269692D0, 0.37220726148105D0, 0.372207010048654D0, 0.372207023598475D0, 0.372207061301007D0, 0.372207102495128D0, 0.372207156735844D0, 0.372207221096294D0, 0.372207297040354D0, 0.372207381990638D0, 0.372207478879362D0, 0.372207588044213D0, 0.372207708503128D0, 0.372207842584591D0, 0.372207991237274D0, 0.372208156448824D0, 0.372208340466725D0, 0.372208544109508D0, 0.37220877068467D0, 0.372209023803806D0, 0.372209305086145D0, 0.372209620875341D0, 0.372209966023574D0, 0.372210367353444D0, 0.372210819041779D0, 0.372211327462498D0, 0.372211903558519D0, 0.372212549102313D0, 0.372213301850905D0, 0.372214119184798D0, 0.372215050211841D0, 0.372216130143523D0, 0.372218084881813D0, 0.372219465012834D0, 0.372221065844759D0/)
avg_earnings_ss_ig(2,2,:)=(/0.554383270087853D0, 0.554410917644114D0, 0.554425915977084D0, 0.5544521327373D0, 0.554496711121797D0, 0.554546985781458D0, 0.554614528015536D0, 0.554671865994041D0, 0.554714205379653D0, 0.554722756976224D0, 0.554730173965784D0, 0.554736994268489D0, 0.554743897530824D0, 0.554743093487189D0, 0.55474578840455D0, 0.55474779305203D0, 0.554746624000803D0, 0.554748796032126D0, 0.554747494449363D0, 0.554747521966766D0, 0.55474931650177D0, 0.554750464693091D0, 0.554751924584304D0, 0.554751509566045D0, 0.5547493319284D0, 0.554749199057982D0, 0.554747548479243D0, 0.554747849814152D0, 0.554747221769678D0, 0.554747793467413D0, 0.554748968061888D0, 0.554750274332522D0, 0.55475094713877D0, 0.554751478774602D0, 0.55475188707038D0, 0.554752886389382D0, 0.55475288352538D0, 0.554754722075539D0, 0.554756266837395D0, 0.554757652374092D0, 0.554758817162375D0, 0.554759563528457D0, 0.554760308228649D0, 0.554760806241801D0, 0.554761693438617D0, 0.554762041318995D0, 0.554762302436727D0, 0.554762658452416D0, 0.554762837231513D0, 0.554762973540029D0, 0.554763081545648D0, 0.554762915253812D0, 0.554762804676444D0, 0.554762627240749D0, 0.554762490200769D0, 0.554762271507659D0, 0.55476207524138D0, 0.55476181817767D0, 0.55476164231792D0, 0.554761462473494D0, 0.554761321963493D0, 0.554761164299291D0, 0.554761020817291D0, 0.554760910259447D0, 0.554760832199799D0, 0.5547608015479D0, 0.554760774682059D0, 0.554760774473158D0, 0.554760797051145D0, 0.554760839439438D0, 0.554760902034709D0, 0.554760979531293D0, 0.554761074201289D0, 0.554761182283857D0, 0.554761304346319D0, 0.554761441869499D0, 0.554761593455363D0, 0.554761760360121D0, 0.554761942513335D0, 0.554762147916329D0, 0.554762374253677D0, 0.554762604419559D0, 0.55476288499004D0, 0.554763198206761D0, 0.554763547575023D0, 0.554763937681398D0, 0.55476437528872D0, 0.554764835304539D0, 0.554765388500174D0, 0.554766014440886D0, 0.554766716324384D0, 0.554767508308589D0, 0.554768407433957D0, 0.554769419648729D0, 0.554770581043158D0, 0.554771928602115D0, 0.554773455178694D0, 0.554775080554106D0, 0.554777182807567D0, 0.554779625294062D0/)
avg_earnings_ss_ig(1,3,:)=(/0.432120687655043D0, 0.432103280825511D0, 0.432119283806918D0, 0.432161439272538D0, 0.43223519178783D0, 0.432310895967989D0, 0.432382243743928D0, 0.432449615497773D0, 0.432502319668361D0, 0.432535754893632D0, 0.432567112053261D0, 0.432593223225261D0, 0.432616053095208D0, 0.432643907372987D0, 0.432663546089506D0, 0.432681992481821D0, 0.432700383069785D0, 0.43271719951733D0, 0.432733956869789D0, 0.432758825951394D0, 0.432774190847745D0, 0.432788962501224D0, 0.432802236086716D0, 0.432815345412941D0, 0.432828435616562D0, 0.432839369049825D0, 0.432850834594581D0, 0.432862390634966D0, 0.432870040914936D0, 0.432877820614077D0, 0.432883449336428D0, 0.432888508929689D0, 0.432893066389995D0, 0.432896673171541D0, 0.432900578450672D0, 0.432903039181644D0, 0.432905835129458D0, 0.432908452698086D0, 0.432910791778727D0, 0.432912577129188D0, 0.432914800420835D0, 0.43291603646629D0, 0.432918299931D0, 0.432919509323986D0, 0.432920469843846D0, 0.432921660712882D0, 0.432922379404584D0, 0.432922844419782D0, 0.432923095044185D0, 0.432923199552284D0, 0.432923263687123D0, 0.432923344080991D0, 0.432923296816492D0, 0.432923172509012D0, 0.432923008348231D0, 0.432922807600814D0, 0.432922614301052D0, 0.432922413809222D0, 0.432922232834199D0, 0.432922066197552D0, 0.43292192218688D0, 0.432921783156547D0, 0.432921658042918D0, 0.432921565296968D0, 0.432921495379284D0, 0.432921452204682D0, 0.432921429842227D0, 0.432921429294262D0, 0.432921442153385D0, 0.432921473244136D0, 0.43292151943101D0, 0.43292158150609D0, 0.432921656366D0, 0.432921744555622D0, 0.432921843099424D0, 0.432921955745083D0, 0.432922080419464D0, 0.432922221309113D0, 0.432922375018134D0, 0.43292254529342D0, 0.432922733706919D0, 0.432922943548229D0, 0.432923175230905D0, 0.432923436022216D0, 0.432923723663686D0, 0.43292404172408D0, 0.432924399701312D0, 0.432924790987957D0, 0.432925244547436D0, 0.432925755116999D0, 0.432926330933289D0, 0.432926982572062D0, 0.432927715438232D0, 0.43292856023026D0, 0.432929492405146D0, 0.432930559648964D0, 0.432931780784143D0, 0.43293543621968D0, 0.432937188429016D0, 0.432939026927261D0/)
avg_earnings_ss_ig(2,3,:)=(/0.634775367104697D0, 0.634775307194339D0, 0.634771288641029D0, 0.634779099468862D0, 0.634813033751434D0, 0.634853725223466D0, 0.634912679619666D0, 0.634958744466581D0, 0.634990843786346D0, 0.634990528733853D0, 0.634984628089209D0, 0.634982587180486D0, 0.634979366777049D0, 0.634975348752059D0, 0.634970185047957D0, 0.634968191709076D0, 0.634961697519229D0, 0.634958900575863D0, 0.63495690459326D0, 0.63495212857744D0, 0.634952171493881D0, 0.634949028571916D0, 0.634944873128483D0, 0.63494391154103D0, 0.634942328981418D0, 0.63494139788128D0, 0.634940033074367D0, 0.634938641676314D0, 0.63493800966674D0, 0.634936552837029D0, 0.634936439762161D0, 0.634937403561113D0, 0.634937360829607D0, 0.634938315888317D0, 0.634938796660522D0, 0.634939011460897D0, 0.634940842811601D0, 0.634941826963557D0, 0.634943162964698D0, 0.634944858650696D0, 0.634945797149635D0, 0.634946516617991D0, 0.634947119151695D0, 0.634947440408136D0, 0.634948056293859D0, 0.634948748443402D0, 0.634949466976894D0, 0.634949964852544D0, 0.634950246290451D0, 0.634950377107271D0, 0.634950472356183D0, 0.634950270817308D0, 0.634950228792779D0, 0.634950131167029D0, 0.634949950354475D0, 0.634949711195714D0, 0.634949460590465D0, 0.634949210400328D0, 0.634948980737471D0, 0.634948770970111D0, 0.634948588147602D0, 0.63494841538291D0, 0.634948266143492D0, 0.634948151667529D0, 0.634948066278919D0, 0.634948018254832D0, 0.63494799886478D0, 0.634948007872555D0, 0.634948037238949D0, 0.634948094493755D0, 0.634948170840665D0, 0.63494826754984D0, 0.634948379209908D0, 0.634948508995932D0, 0.634948651013795D0, 0.634948808864468D0, 0.634948983154312D0, 0.634949175658168D0, 0.634949387499099D0, 0.634949621460165D0, 0.634949878185281D0, 0.634950162984658D0, 0.634950477468732D0, 0.634950827014147D0, 0.634951217259417D0, 0.634951652937534D0, 0.634952141295173D0, 0.634952689767343D0, 0.634953303977568D0, 0.634953995726919D0, 0.634954774154834D0, 0.634955650943266D0, 0.634956644526325D0, 0.63495777565245D0, 0.634959058177416D0, 0.634960538655708D0, 0.634962246792513D0, 0.634964094333514D0, 0.634966401291807D0, 0.634969059564426D0/)
avg_earnings_ss_ig(1,4,:)=(/0.503124304312219D0, 0.503073285748494D0, 0.503072228390574D0, 0.503100901491108D0, 0.503157573220025D0, 0.503216509411065D0, 0.503276126439433D0, 0.503333477645308D0, 0.503378691967356D0, 0.503399815087156D0, 0.503424558311247D0, 0.503443926320858D0, 0.503462837385526D0, 0.503479440881293D0, 0.503495577087707D0, 0.503512296882627D0, 0.503525493008215D0, 0.503537180613618D0, 0.503550690207083D0, 0.503565115954079D0, 0.503575028401225D0, 0.503588544377982D0, 0.503598177355135D0, 0.503611454383143D0, 0.503622297859845D0, 0.503634271870182D0, 0.503646256317376D0, 0.503654377784601D0, 0.503662897286852D0, 0.50366914114028D0, 0.503674390284371D0, 0.503679072706608D0, 0.50368242171022D0, 0.50368664407263D0, 0.503691871189968D0, 0.503695266608694D0, 0.503698062627068D0, 0.503699967650262D0, 0.503701379410164D0, 0.503703051648D0, 0.503705336392465D0, 0.503707856731665D0, 0.503709212658341D0, 0.503710541596315D0, 0.503712376809902D0, 0.503713247232161D0, 0.503714063965732D0, 0.503714753558287D0, 0.503715164195431D0, 0.503715413820802D0, 0.503715849732228D0, 0.503716051133883D0, 0.503716003851384D0, 0.503715855902679D0, 0.503715669704793D0, 0.503715450396741D0, 0.503715219556781D0, 0.503714994042647D0, 0.503714786015003D0, 0.503714595622033D0, 0.503714425331639D0, 0.503714262394721D0, 0.503714116386518D0, 0.503714002906014D0, 0.503713917400934D0, 0.503713863636683D0, 0.503713833337741D0, 0.503713828117204D0, 0.503713840086538D0, 0.503713874438299D0, 0.503713926997084D0, 0.503713998684461D0, 0.503714084145965D0, 0.503714185910543D0, 0.50371429967711D0, 0.50371442866726D0, 0.503714571356157D0, 0.503714728851166D0, 0.503714903281028D0, 0.503715096254495D0, 0.503715308759987D0, 0.503715544654598D0, 0.50371580513737D0, 0.503716095696106D0, 0.503716417404365D0, 0.503716778458662D0, 0.503717181433796D0, 0.503717624189486D0, 0.5037181352377D0, 0.50371871242924D0, 0.503719364036055D0, 0.503720098155353D0, 0.503720927887715D0, 0.503721887876793D0, 0.503722950889588D0, 0.503724169619512D0, 0.503725559148607D0, 0.503728865073383D0, 0.503730797842659D0, 0.503732889128298D0/)
avg_earnings_ss_ig(2,4,:)=(/0.719054387108188D0, 0.719041797178325D0, 0.719017873212277D0, 0.719026093862415D0, 0.719057147626552D0, 0.71909132483494D0, 0.719141711434138D0, 0.719189398546082D0, 0.719222813305749D0, 0.719216736227835D0, 0.71920832211704D0, 0.719197812649659D0, 0.719189522089343D0, 0.719180265549008D0, 0.719175509074994D0, 0.719172142634119D0, 0.71916898760434D0, 0.719165298943024D0, 0.719160736382437D0, 0.719156240793699D0, 0.719154975145529D0, 0.719150661374486D0, 0.719146976639779D0, 0.719144152728936D0, 0.719141594791261D0, 0.719139445724539D0, 0.719136651305541D0, 0.719134679509652D0, 0.719134403414069D0, 0.719134127210941D0, 0.719134756939893D0, 0.719134000259197D0, 0.719133653686856D0, 0.719133850917486D0, 0.719135385619719D0, 0.719136417821114D0, 0.719138091596901D0, 0.719138515672924D0, 0.719139295423155D0, 0.719140444633616D0, 0.719141920164707D0, 0.719142455270087D0, 0.719142754942445D0, 0.719143568841237D0, 0.719144125315476D0, 0.719144829196443D0, 0.719145475507875D0, 0.719145984662304D0, 0.719146329063705D0, 0.719146548267627D0, 0.719146582162223D0, 0.719146718186987D0, 0.719146716492382D0, 0.719146589470564D0, 0.71914640441196D0, 0.719146155079076D0, 0.719145887668564D0, 0.71914559983451D0, 0.719145322456227D0, 0.719145100073291D0, 0.719144904855735D0, 0.719144713769486D0, 0.719144593216105D0, 0.719144475209464D0, 0.719144378989029D0, 0.719144320338292D0, 0.719144292165549D0, 0.719144291003887D0, 0.719144319035461D0, 0.719144380523075D0, 0.719144462512948D0, 0.719144568017079D0, 0.71914468895008D0, 0.719144832271216D0, 0.719144990175539D0, 0.71914516514091D0, 0.719145359829179D0, 0.719145576450265D0, 0.719145812324018D0, 0.719146072330132D0, 0.719146356934256D0, 0.719146672631853D0, 0.719147021689185D0, 0.719147409499937D0, 0.719147840614168D0, 0.719148320607945D0, 0.719148858329699D0, 0.719149484364612D0, 0.719150159805419D0, 0.719150919721648D0, 0.719151772632399D0, 0.719152732947275D0, 0.719153816620836D0, 0.719155006006229D0, 0.719156445032321D0, 0.719158041689256D0, 0.719159898975244D0, 0.719161982031455D0, 0.719164489329798D0, 0.719167362028369D0/)
avg_earnings_ss_ig(1,5,:)=(/0.592937917248836D0, 0.592869476497215D0, 0.592835482610043D0, 0.592829597333379D0, 0.59287483425228D0, 0.592918429060269D0, 0.592968479955481D0, 0.593009961384367D0, 0.593044431604541D0, 0.593054713901585D0, 0.593070024433374D0, 0.593075080510352D0, 0.593086503666883D0, 0.593094724485994D0, 0.593101255683551D0, 0.593112335349059D0, 0.593120249906221D0, 0.593128822946254D0, 0.593136653331486D0, 0.593147528954679D0, 0.593157619587468D0, 0.593162918373674D0, 0.593172841450964D0, 0.593183548481601D0, 0.593190810334912D0, 0.593196520025753D0, 0.593206402747825D0, 0.593215338509014D0, 0.593222822974669D0, 0.593227524981121D0, 0.593232807354101D0, 0.593237969457878D0, 0.593239460615596D0, 0.593242063515078D0, 0.593243070416581D0, 0.593246005448441D0, 0.593248318895707D0, 0.593250980784265D0, 0.593254353617826D0, 0.593255959892491D0, 0.593258058096081D0, 0.593262054453562D0, 0.593263702046279D0, 0.593266107633081D0, 0.59326757582443D0, 0.59326756734062D0, 0.593268511139431D0, 0.593268986224712D0, 0.593269154575622D0, 0.593269383566896D0, 0.593269414090782D0, 0.59326920505527D0, 0.593269014602499D0, 0.593268826071714D0, 0.593268614502748D0, 0.593268386737725D0, 0.593268147710585D0, 0.593267908812968D0, 0.593267690321956D0, 0.593267494330252D0, 0.593267323329713D0, 0.593267160931874D0, 0.59326701855599D0, 0.593266910878878D0, 0.593266834090121D0, 0.593266796357841D0, 0.593266783650637D0, 0.593266797127423D0, 0.593266833207676D0, 0.593266894502909D0, 0.593266973410615D0, 0.593267071395086D0, 0.593267182680525D0, 0.593267310928599D0, 0.593267451951169D0, 0.593267608068872D0, 0.593267779952174D0, 0.593267966594725D0, 0.593268170341374D0, 0.593268392855417D0, 0.593268635301121D0, 0.593268903986743D0, 0.593269200049704D0, 0.593269531902036D0, 0.593269898234833D0, 0.593270307079036D0, 0.593270766650973D0, 0.593271271790852D0, 0.593271848451865D0, 0.593272496538373D0, 0.593273223118131D0, 0.593274043192768D0, 0.593274970431866D0, 0.593276033952418D0, 0.593277219399114D0, 0.593278632693556D0, 0.593280176728673D0, 0.593281832604603D0, 0.59328382474417D0, 0.593286167491816D0/)
avg_earnings_ss_ig(2,5,:)=(/0.79596158165423D0, 0.795924574285371D0, 0.795897433319792D0, 0.79589147203702D0, 0.795914189729399D0, 0.795952082045121D0, 0.796014811703934D0, 0.796056646660726D0, 0.796090847072318D0, 0.796084753845028D0, 0.796073682343483D0, 0.796062734661835D0, 0.796050994368959D0, 0.79604005354336D0, 0.796036611239055D0, 0.796033242583683D0, 0.796025306898574D0, 0.796019195126248D0, 0.79601332816449D0, 0.79600622500725D0, 0.796009217792161D0, 0.796008691771363D0, 0.796009678556434D0, 0.796007386801201D0, 0.796009429788038D0, 0.796008227200914D0, 0.79600861694233D0, 0.796009885779731D0, 0.796010842019294D0, 0.796014013939581D0, 0.796014140067607D0, 0.796016776564859D0, 0.796018722562102D0, 0.796020356669255D0, 0.79602304118798D0, 0.796024970532813D0, 0.796026884310942D0, 0.796027380619678D0, 0.796028329998824D0, 0.79602911695812D0, 0.796030291588826D0, 0.796032522972813D0, 0.796033734872873D0, 0.796033338340882D0, 0.796033840990187D0, 0.796034162229818D0, 0.796034722250225D0, 0.796035161237962D0, 0.796035435915662D0, 0.796035515278542D0, 0.796035492194994D0, 0.796035393896042D0, 0.796035210067784D0, 0.796034988145276D0, 0.796034768025196D0, 0.796034501713051D0, 0.7960342279463D0, 0.796033940543538D0, 0.796033678723741D0, 0.796033446500943D0, 0.796033249506165D0, 0.796033055552381D0, 0.796032920488742D0, 0.796032809740262D0, 0.796032722822472D0, 0.796032679186447D0, 0.796032667561406D0, 0.796032685071754D0, 0.796032737230599D0, 0.796032818946773D0, 0.79603292105932D0, 0.796033047528305D0, 0.796033191852046D0, 0.796033356313315D0, 0.796033540599849D0, 0.796033741545046D0, 0.796033962401586D0, 0.7960342049876D0, 0.796034466554668D0, 0.796034756432256D0, 0.796035073469493D0, 0.796035418700443D0, 0.796035802393578D0, 0.796036221618688D0, 0.79603667974078D0, 0.796037204217376D0, 0.7960377814689D0, 0.796038435037979D0, 0.796039164410348D0, 0.79603997158827D0, 0.79604087178319D0, 0.796041885114801D0, 0.796043028512375D0, 0.796044302554159D0, 0.796045816539323D0, 0.796047522124019D0, 0.796049507026864D0, 0.796051803498904D0, 0.796054473917819D0, 0.796057512764446D0/)




if(routine==4 .or. routine ==2) then
    lambda0=aggregates_end(13)
!    lambda0=.714327013641034D0
end if
!do it=1,typesN
    ss_ig=ss_start*(Ce_for_share_start*(pe_start+taue)+C_for_share_start)/(Ce_for_share_start*pe_start+C_for_share_start)
!end do
if (optimal_baseline==0) then
 OPEN(227, FILE="/<filepath starting steady state>/sim_labor_output.dat",FORM="formatted")
    read(227,887) (sim_labor_start)
    close(227)

 OPEN(227, FILE="/<filepath starting steady state>/sim_energy_output.dat",FORM="formatted")
    read(227,887) (sim_energy_start)
    close(227)

 OPEN(227, FILE="/<filepath starting steady state>/sim_consumption_output.dat",FORM="formatted")
    read(227,887) (sim_consumption_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steady state>/sim_capital_output.dat",FORM="formatted")
    read(227,887) (sim_capital_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steady state>/sim_labor_earnings_output.dat",FORM="formatted")
    read(227,887) (sim_labor_earnings_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steady state>/sim_labor_tax_output.dat",FORM="formatted")
    read(227,887) (sim_labor_tax_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steady state>/sim_capital_earnings_output.dat",FORM="formatted")
    read(227,887) (sim_capital_earnings_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steady state>/sim_capital_tax_output.dat",FORM="formatted")
    read(227,887) (sim_capital_tax_start)
    close(227)

sim_total_after_tax_earnings_start=sim_labor_earnings_start+sim_capital_earnings_start-sim_capital_tax_start-sim_labor_tax_start
else
 OPEN(227, FILE="/<filepath starting steady state>/unconstrained_optimal/sim_labor_output.dat",FORM="formatted")
    read(227,887) (sim_labor_start)
    close(227)

 OPEN(227, FILE="/<filepath starting steady state>/unconstrained_optimal/sim_energy_output.dat",FORM="formatted")
    read(227,887) (sim_energy_start)
    close(227)

 OPEN(227, FILE="/<filepath starting steady state>/unconstrained_optimal/sim_consumption_output.dat",FORM="formatted")
    read(227,887) (sim_consumption_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steady state>/unconstrained_optimal/sim_capital_output.dat",FORM="formatted")
    read(227,887) (sim_capital_start)
    close(227)

 OPEN(227, FILE="/<filepath starting steady state>/unconstrained_optimal/sim_labor_earnings_output.dat",FORM="formatted")
    read(227,887) (sim_labor_earnings_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steady state>/unconstrained_optimal/sim_labor_tax_output.dat",FORM="formatted")
    read(227,887) (sim_labor_tax_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steady state>/unconstrained_optimal/sim_capital_earnings_output.dat",FORM="formatted")
    read(227,887) (sim_capital_earnings_start)
    close(227)
 OPEN(227, FILE="/<filepath starting steady state>/unconstrained_optimal/sim_capital_tax_output.dat",FORM="formatted")
    read(227,887) (sim_capital_tax_start)
    close(227)

sim_total_after_tax_earnings_start=sim_labor_earnings_start+sim_capital_earnings_start-sim_capital_tax_start-sim_labor_tax_start
end if


!change
if(routine==2 .or. routine == 3 .or. routine ==4 .or. routine==5) then
    lump_sum_ig=0.00D0
else if(routine==1) then
    lump_sum_ig=0.05D0
else
    print*,"Not Specified Right"
end if



!Set Up population for simulations
seed=53415691
do s=1,simulationsN
do i=1,J
    random(i,s)=randu(seed)
end do
end do

i=1
do s=1,typesN
    sim_types(1,i:i+int(simulationsN/(typesN))-1,1)=s
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(typesN))-1
end if
    i=i+int(simulationsN/(typesN))
end do

!sim_types(1,1:int((simulationsN)/2))=1
!sim_types(1,int((simulationsN)/2)+1:simulationsN)=2

do s=2,J
    sim_types(s,:,1)=sim_types(1,:,1)
end do


seed=52345691
do s=1,simulationsN
do i=1,J
    random(i,s)=randu(seed)
end do
end do


i=1
do it = 1,typesN
    sim_shocks(1,i:i+int(simulationsN/(2*(typesN)))-1,1)=shocks_start1
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(2*(typesN)))-1
end if
    i=i+int(simulationsN/(2*(typesN)))
    sim_shocks(1,i:i+int(simulationsN/(2*(typesN)))-1,1)=shocks_start2
if (my_id==1) then
print*,"i i+",i,i+int(simulationsN/(2*(typesN)))-1
end if
    i=i+int(simulationsN/(2*(typesN)))
end do

do i=1,shocksN
do s=1,shocksN
do it1=1,Jnr-1
do it2=1,2
    shock_cumpi(i,s,it1,it2)=sum(shock_pi(i,1:s,it1,it2))
end do
end do
end do
end do

do it1=Jnr,J
do it2=1,2
    shock_cumpi(:,:,it1,it2)=shock_cumpi(:,:,Jnr-1,it2)
end do
end do

shock_cumpi(:,shocksN,:,:)=1.0D0


!Low
do s=1,int(simulationsN/2)
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s,1)=minloc(shock_cumpi(sim_shocks(i-1,s,1),:,i,1),1,mask=shock_cumpi(sim_shocks(i-1,s,1),:,i,1).ge.random(i,s))
    else
        sim_shocks(i,s,1)=sim_shocks(i-1,s,1)
    end if
end do
end do
!High
do s=int(simulationsN/2)+1,simulationsN
do i =2,J
    if (i<Jnr) then
        sim_shocks(i,s,1)=minloc(shock_cumpi(sim_shocks(i-1,s,1),:,i,2),1,mask=shock_cumpi(sim_shocks(i-1,s,1),:,i,2).ge.random(i,s))
    else
        sim_shocks(i,s,1)=sim_shocks(i-1,s,1)
    end if
end do
end do

if(my_id .eq. root_process) then
print*,"I got herea"
end if
ss_count=0.0D0
do it =1 ,simulationsN
    if(sim_shocks(Jnr,it,1).le. shocksN/2) then
        ss_count(sim_types(1,it,1),sim_shocks(Jnr,it,1),1)=ss_count(sim_types(1,it,1),sim_shocks(Jnr,it,1),1)+1.0D0
    else
        ss_count(sim_types(1,it,1),sim_shocks(Jnr,it,1)-shocksN/2,1)=ss_count(sim_types(1,it,1),sim_shocks(Jnr,it,1)-shocksN/2,1)+1.0D0
    end if
end do
if(my_id .eq. root_process) then
print*,"I got hereb"
end if

do it=2,yr_T
    ss_count(:,:,it)=ss_count(:,:,1)
    sim_shocks(:,:,it)=sim_shocks(:,:,1)
    sim_types(:,:,it)=sim_types(:,:,1)
end do

!MASTER
if (my_id .eq. root_process) then


if(my_id .eq. root_process) then
print*,"I got here"
end if



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ENTER FUNCTION HERE IF MAKING A FUNCTION
!Initial guesses
K=K_ig
N=N_ig
tr=tr_ig
Ec=Ec_ig
do it=1,yr_T
    ss(:,:,it)=ss_ig
end do
lump_sum=lump_sum_ig
avg_earnings_ss=avg_earnings_ss_ig

if(my_id .eq. root_process) then
print*,"here3",rebate_level,rebate_slope,lambda1,K
end if

!!!!!!!!!!!!!!!!!!!!!!!! Initialize loop stuff
tax_loop = 1
tax_diff=1.0D0
tax_diff_one=1.0D0
non_clear=0
bad_combination_inner_stop=0
bad_combination_inner_nostop=0
bad_combination_outer=0
diff_parameter=1.0D0
out_loop=1
allocate(zeros2(max_outer_loops+1))
zeros2=0.0D0
do it=1,yr_T
    taxes_clear(:,it)=zeros2
end do
deallocate(zeros2)
bad_combination_inner_stop=0
bad_combination_outer=0

!Loop for price convergence
DO WHILE (diff_parameter>tol)
   tax_loop = 1
   tax_diff_one=1.0D0
   non_clear=0
DO WHILE (tax_diff_one >tol)

    do yr=yr_T,1,-1
    year_calculating=yr
    print*,"yr",yr,out_loop

    passed_parameters(1)=K(yr)
    passed_parameters(2)=N(yr)
    passed_parameters(3)=Ec(yr)
    passed_parameters(4)=Atfp
    passed_parameters(5)=Ae
    passed_parameters(6)=taue

print*,"parameters_passed",passed_parameters

        allocate(share_guesses(6))
        allocate(share_guesses2(8))
        allocate(share_guesses3(8))
        if(out_loop==1 ) then
            share_guesses=(/.025D0,.035D0,.05D0,.015D0,.04D0,.01D0/)
        else
            share_guesses=(/.5D0,Ep(yr),Ep(yr)*1.1D0,Ep(yr)*.8D0,Ep(yr)*1.05D0,Ep(yr)*.9D0/)
        end if
        share_guesses2=(/.5D0,.65D0,.8D0,.9D0,.94D0,.95D0,.96D0,.98D0/)
        share_guesses3=(/.5D0,.65D0,.8D0,.9D0,.94D0,.95D0,.96D0,.98D0/)
        do share_it=1,6
        do share_it2=1,8
        do share_it3=1,8
            xguesses_check(1)=share_guesses(share_it)
            xguesses_check(2)=K(yr)*share_guesses2(share_it2)
            xguesses_check(3)=N(yr)*share_guesses3(share_it3)
            CALL energy_root_finder(xguesses_check,check_dummy_out,3)
            if(share_it==1 .and. share_it2==1 .and. share_it3==1) then
                check_dummy=check_dummy_out(1)**2.0D0+check_dummy_out(2)**2.0D0+check_dummy_out(3)**2.0D0
                xguess(1)=xguesses_check(1)
                xguess(2)=xguesses_check(2)
                xguess(3)=xguesses_check(3)
            elseif (check_dummy>(check_dummy_out(1)**2.0D0+check_dummy_out(2)**2.0D0+check_dummy_out(3)**2.0D0)) then
                xguess(1)=xguesses_check(1)
                xguess(2)=xguesses_check(2)
                xguess(3)=xguesses_check(3)
                check_dummy=check_dummy_out(1)**2.0D0+check_dummy_out(2)**2.0D0+check_dummy_out(3)**2.0D0
            end if
        end do
        end do
        end do
        deallocate(share_guesses)
        deallocate(share_guesses2)
        deallocate(share_guesses3)
print*,"xguess in",xguess

    CALL D_NEQNF(energy_root_finder,solver_out,ERRREL=0.0000001D0,xguess=xguess)
print*,"xguess out",solver_out

    Ep(yr)=solver_out(1)
    E(yr)=Ep(yr)+Ec(yr)
    K1(yr)=solver_out(2)
    K2(yr)=K(yr)-K1(yr)
    N1(yr)=solver_out(3)
    N2(yr)=N(yr)-N1(yr)


    Y(yr)=Atfp*(K1(yr)**(alpha1)*N1(yr)**(1.0D0-alpha1-theta))*Ep(yr)**theta
    pe(yr)=Atfp*theta*(K1(yr)/Ep(yr))**alpha1*(N1(yr)/Ep(yr))**(1.0D0-alpha1-theta)-taue
    w(yr)=Atfp*(1.0D0-alpha1-theta)*(K1(yr)/N1(yr))**(alpha1)*(Ep(yr)/N1(yr))**theta
    r(yr)=Atfp*alpha1*(N1(yr)/K1(yr))**(1.0D0-alpha1-theta)*(Ep(yr)/K1(yr))**theta-delta
    tauss(yr)=sum(sum(ss(:,:,yr)*ss_count(:,:,yr),2),1)/real(simulationsN)*sum(population(Jnr:J,2))/&
        (sum(sum(avg_earnings_ss(:,:,yr)*ss_count(:,:,yr),2),1)/real(simulationsN)*sum(population(1:Jnr-1,2)))

    ss_cap(yr)=avg_earnings(yr)*2.42
print*,"prices",yr,w(yr),r(yr),avg_earnings_scale,pe(yr),tauss(yr),lambda0(yr),lambdak(yr)

!Extra loop for convergence on tax parameters (if turned on need tax parameter convergence begore prices update)
!Solve for policy funcation
        v(:,:,:,:,yr) = -1.0D0/0.0D0
        v(:,:,:,J+1-jretired,yr)=0.0D0
        v_retired(:,:,:,:,yr) = -1.0D0/0.0D0
        kprime_dummy=0
        v_retired(:,:,capital_positive:kgridN,J+1-jretired,yr)=0.0D0

!Last generation
            do it=capital_positive,kgridN
                do it2=1,typesN
                    do e_it=1,shocksN
                        if(e_it .le. (shocksN)/2) then
                            ctot80(it,1)=retcon(kgrid(it),0.0D0,tr(yr),r(yr),ss(it2,e_it,yr),avg_earnings(yr),lambdak(yr))
                            energy_miss=energy_sharer(0.00D0,ctot80(it,1)*.15D0,ctot80(it,1),energy_share_equation,.0001D0,c80(it,1),ctot80(it,1),J)
                            e80(it,1)=(ctot80(it,1)-c80(it,1))/(pe(yr)+taue)
                            v_retired(it2,e_it,it,J-jretired,yr) =utility_last(c80(it,1),e80(it,1),adult_equivalence(J))


                        else
                            v_retired(it2,e_it,it,J-jretired,yr) =v_retired(it2,e_it-(shocksN/2),it,J-jretired,yr)
                        end if
                    end do
                end do
            end do


!Retired generations
        allocate(zeros(kgridN,1))
        zeros=0.0D0
        do s =(J-1),Jnr, -1
            do e_it=1,shocksN
            if(e_it .le. (shocksN)/2) then
                do it=1,typesN
                do kit=capital_positive, kgridN
                    con=-1.0D0
                    con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,kgrid(kit),kgrid(capital_positive:),&
                        tr(yr),r(yr),ss(it,e_it,yr),avg_earnings(yr),lambdak(yr))
                    con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))
                    if (con(kgridN)>0.015+ebar) then
                        con_pos=kgridN
                    end if
                    if (con(capital_positive+1)<=0.015+ebar) then
                        v_retired(it,e_it,kit,s-jretired,yr)=-1.0D0/0.0D0
                    else
                    CALL retired_valuer(kgrid(kit),0.0D0,kgrid(con_pos),ss(it,e_it,yr),tr(yr),r(yr),survive(s),v_retired(it,e_it,:,s+1-jretired,yr+1),&
                        v_retired(it,e_it,kit,s-jretired,yr),kprime_dummy,avg_earnings(yr),s,lambdak(yr))
                    end if
               end do
               end do
            else
                do it=1,typesN
                do kit=capital_positive, kgridN
                    v_retired(it,e_it,kit,s-jretired,yr)=v_retired(it,e_it-(shocksN/2),kit,s-jretired,yr)
                end do
                end do
            end if
            end do
        end do
        deallocate(zeros)

if(my_id==0) then
print*,"here4"
end if


!Working generations
v(:,:,:,Jnr,yr+1)=v_retired(:,:,:,Jnr-jretired,yr+1)

                prices(1)=tr(yr)
                prices(2)=lambda0(yr)
                prices(3)=w(yr)
                prices(4)=r(yr)
                prices(5)=tauss(yr)
                prices(6)=rebate_level(yr)
                prices(7)=ss_cap(yr)
                prices(8)=lump_sum(yr)
                prices(9)=avg_earnings(yr)
                prices(10)=lambdak(yr)
                prices(11)=pe(yr)
                CALL MPI_Bcast( prices, 11, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
                CALL MPI_Bcast( ss(:,:,yr), typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
        do s =(Jnr-1),2,-1
            CALL MPI_Bcast( v(:,:,:,s+1,yr+1),typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            end_row=0
            do an_id = 1, num_procs_run -1
                start_row=end_row+1
                if (an_id>large_inc) then
                    end_row=start_row+increment-1
                else
                    end_row=start_row+increment
                end if
                call MPI_SEND(start_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
                call MPI_SEND(end_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
            end do
            do an_id=1,num_procs_run -1
                call MPI_RECV( start_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( end_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
                call MPI_RECV( v(:,:,start_row:end_row,s,yr),typesN*shocksN*(end_row-start_row+1) , MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            end do

        end do
        s=1
        CALL MPI_Bcast( v(:,:,:,s+1,yr+1),typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )


if(my_id==0) then
print*,"here5"
end if

!yr loop
end do

if(my_id==0) then
print*,"value function done"
end if
!! Simulate peoples

    CALL MPI_Bcast( v_retired(:,:,:,:,:),typesN*shocksN*kgridN*(J+1-jretired)*(yr_T+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end_row=0
    do an_id = 1, num_procs_run2 -1
        start_row=end_row+1
        if (an_id>large_inc2) then
            end_row=start_row+increment2-1
        else
            end_row=start_row+increment2
        end if
        call MPI_SEND(start_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
        call MPI_SEND(end_row, 1, MPI_INTEGER, an_id, my_id, MPI_COMM_WORLD, ierr)
    end do

    do an_id=1,num_procs_run2 -1
        call MPI_RECV( start_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( end_row,1 , MPI_INTEGER, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_labor(:,start_row:end_row,:),(J)*(end_row-start_row+1)*(yr_T), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_earnings(:,start_row:end_row,:),(J)*(end_row-start_row+1)*(yr_T), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_capital(:,start_row:end_row,:),(J+1)*(end_row-start_row+1)*(yr_T+1), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_consumption(:,start_row:end_row,:),(J)*(end_row-start_row+1)*(yr_T), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_energy(:,start_row:end_row,:),(J)*(end_row-start_row+1)*(yr_T), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_v(start_row:end_row,:),(end_row-start_row+1)*(yr_T), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_total_income(:,start_row:end_row,:),(J)*(end_row-start_row+1)*(yr_T), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
        call MPI_RECV( sim_after_tax_income(:,start_row:end_row,:),(J)*(end_row-start_row+1)*(yr_T), MPI_DOUBLE_PRECISION, an_id, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

    end do


do yr=1,yr_T
year_calculating=yr
!Calculate Aggregate Profiles

agg_labor(:,1,yr)=sum(sim_labor(:,1:int((simulationsN)/2),yr),2)/int((simulationsN)/2.0D0)
agg_labor(:,2,yr)=sum(sim_labor(:,int((simulationsN)/2)+1:simulationsN,yr),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_capital(:,1,yr)=sum(sim_capital(:,1:int((simulationsN)/2),yr),2)/int((simulationsN)/2.0D0)
agg_capital(:,2,yr)=sum(sim_capital(:,int((simulationsN)/2)+1:simulationsN,yr),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_consumption(:,1,yr)=sum(sim_consumption(:,1:int((simulationsN)/2),yr),2)/int((simulationsN)/2.0D0)
agg_consumption(:,2,yr)=sum(sim_consumption(:,int((simulationsN)/2)+1:simulationsN,yr),2)/(simulationsN-int((simulationsN)/2.0D0))
agg_energy(:,1,yr)=sum(sim_energy(:,1:int((simulationsN)/2),yr),2)/int((simulationsN)/2.0D0)
agg_energy(:,2,yr)=sum(sim_energy(:,int((simulationsN)/2)+1:simulationsN,yr),2)/(simulationsN-int((simulationsN)/2.0D0))

if(yr<5) then
print*,"year",yr,agg_labor(:,1,yr),agg_capital(:,1,yr),agg_consumption(:,1,yr)
print*,"lab,cap,con",sim_labor(:,1,yr),sim_capital(:,1,yr),sim_consumption(:,1,yr)
end if

ind_eng_tax(yr)=0.0D0
do an_id=1,J
    ind_eng_tax(yr)=ind_eng_tax(yr)+taue*(agg_energy(an_id,1,yr)*.5D0+.5D0*agg_energy(an_id,2,yr))*population(an_id,2)
end do
!!Find lambda 0 that clears the market


!1 lump_sum_rebate
!2 throw_away
!3 cap_tax_rebate
!4 labor_rebate

if (routine == 3) then
CALL lambdak_optimizer(lambda0_1(yr),gov_spend(yr)-Ep(yr)*taue,sim_labor(:,:,yr),sim_energy(:,:,yr),&
    sim_shocks(:,:,yr),sim_types(:,:,yr),sim_capital(:,:,yr),hbar,w(yr),shocks,&
    types,partprod,r(yr),tr(yr),population(:,2),avg_earnings(yr),ss(:,:,yr),survive,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))




elseif (routine ==4) then
!offset other taxes
CALL lambda0_optimizer(lambda0_1(yr),gov_spend(yr)-Ep(yr)*taue,sim_labor(:,:,yr),sim_energy(:,:,yr),&
    sim_shocks(:,:,yr),sim_types(:,:,yr),sim_capital(:,:,yr),hbar,w(yr),shocks,&
    types,partprod,r(yr),tr(yr),population(:,2),avg_earnings(yr),ss(:,:,yr),survive,ss_cap(yr),lambda1(yr),lambdak(yr))

elseif (routine ==5) then
!offset other taxes
CALL rebate_optimizer(rebate_level1(yr),gov_spend(yr)-Ep(yr)*taue,sim_labor(:,:,yr),sim_energy(:,:,yr),&
    sim_shocks(:,:,yr),sim_types(:,:,yr),sim_capital(:,:,yr),hbar,w(yr),shocks,&
    types,partprod,r(yr),tr(yr),population(:,2),avg_earnings(yr),ss(:,:,yr),survive,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))!rebate



else if( routine==1  ) then
allocate(dummy_zeros(J,simulationsN))
dummy_zeros=0.0D0
CALL lambda0_optimizer(lambda0_1(yr),gov_spend(yr),sim_labor(:,:,yr),dummy_zeros,sim_shocks(:,:,yr),sim_types(:,:,yr),&
    sim_capital(:,:,yr),hbar,w(yr),shocks,&
    types,partprod,r(yr),tr(yr),population(:,2),avg_earnings(yr),ss(:,:,yr),survive,ss_cap(yr),lambda1(yr),lambdak(yr))
else if(routine==2) then
    lambda0_1(yr)=lambda0(yr)
else
    print*,"Not Specified Right"
end if

        if (routine == 3) then
            tax_diff(yr)=abs((lambda0_1(yr))-lambdak(yr))/lambdak(yr)
            tax_diff_keep(yr)=tax_diff(yr)
print*,"tax diff",tax_diff(yr),lambda0_1(yr),lambdak(yr)
            lambdak(yr)=(a1*lambda0_1(yr)+(1.0D0-a1)*lambdak(yr))
        elseif(routine==5) then
            tax_diff(yr)=abs((rebate_level1(yr))-rebate_level(yr))/rebate_level(yr)
            tax_diff_keep(yr)=tax_diff(yr)
print*,"tax diff",tax_diff(yr),rebate_level1(yr),rebate_level(yr)
            rebate_level(yr)=(a1*rebate_level1(yr)+(1.0D0-a1)*rebate_level(yr))
        elseif(routine==4) then
        govt_taxes1(yr)=tax_differ_lambda0(lambda0(yr),sim_labor(:,:,yr),sim_energy(:,:,yr),sim_shocks(:,:,yr),sim_types(:,:,yr),&
            sim_capital(:,:,yr),hbar,w(yr),shocks,types,partprod,&
            r(yr),tr(yr),gov_spend(yr)-Ep(yr)*taue,population(:,2),avg_earnings(yr),ss(:,:,yr),survive,&
            ss_cap(yr),lambda1(yr),lambdak(yr))
        govt_taxes2(yr)=tax_differ_lambda0(lambda0_1(yr),sim_labor(:,:,yr),sim_energy(:,:,yr),sim_shocks(:,:,yr),sim_types(:,:,yr),&
            sim_capital(:,:,yr),hbar,w(yr),shocks,types,partprod,&
            r(yr),tr(yr),gov_spend(yr)-Ep(yr)*taue,population(:,2),avg_earnings(yr),ss(:,:,yr),survive,&
            ss_cap(yr),lambda1(yr),lambdak(yr))
        govt_taxes3(yr)=tax_differ_lambda0(0.4D0,sim_labor(:,:,yr),sim_energy(:,:,yr),sim_shocks(:,:,yr),sim_types(:,:,yr),&
            sim_capital(:,:,yr),hbar,w(yr),shocks,types,partprod,&
            r(yr),tr(yr),gov_spend(yr)-Ep(yr)*taue,population(:,2),avg_earnings(yr),ss(:,:,yr),survive,&
            ss_cap(yr),lambda1(yr),lambdak(yr))
        govt_taxes4(yr)=tax_differ_lambda0(.90D0,sim_labor(:,:,yr),sim_energy(:,:,yr),sim_shocks(:,:,yr),sim_types(:,:,yr),&
            sim_capital(:,:,yr),hbar,w(yr),shocks,types,partprod,&
            r(yr),tr(yr),gov_spend(yr)-Ep(yr)*taue,population(:,2),avg_earnings(yr),ss(:,:,yr),survive,&
            ss_cap(yr),lambda1(yr),lambdak(yr))
!            end if

            tax_diff(yr)=abs((lambda0_1(yr))-lambda0(yr))/lambda0(yr)
            tax_diff_keep(yr)=tax_diff(yr)
print*,"tax diff",yr,tax_diff(yr),lambda0_1(yr),lambda0(yr),gov_spend(yr)-Ep(yr)*taue,govt_taxes1(yr),govt_taxes2(yr),govt_taxes3(yr),govt_taxes4(yr)
            lambda0(yr)=(a1*lambda0_1(yr)+(1.0D0-a1)*lambda0(yr))
        elseif(routine==2) then
            tax_diff(yr)=0.0D0
            tax_diff_keep(yr)=tax_diff(yr)
        else
            tax_diff(yr)=abs((lambda0_1(yr))-lambda0(yr))/lambda0(yr)
            tax_diff_keep(yr)=tax_diff(yr)
print*,"tax diff",yr,tax_diff(yr),lambda0_1(yr),lambda0(yr),gov_spend(yr)-Ep(yr)*taue
            lambda0(yr)=(a1*lambda0_1(yr)+(1.0D0-a1)*lambda0(yr))
        end if

if (routine == 3 .or. routine ==4 ) then
!offset other taxes
        govt_taxes(yr)=tax_differ_lambda0(lambda0(yr),sim_labor(:,:,yr),sim_energy(:,:,yr),sim_shocks(:,:,yr),sim_types(:,:,yr),&
            sim_capital(:,:,yr),hbar,w(yr),shocks,types,partprod,&
            r(yr),tr(yr),gov_spend(yr)-Ep(yr)*taue,population(:,2),avg_earnings(yr),ss(:,:,yr),survive,&
            ss_cap(yr),lambda1(yr),lambdak(yr))
else if (routine==5) then
!offset other taxes
        govt_taxes(yr)=tax_differ_rebate(rebate_level(yr),sim_labor(:,:,yr),sim_energy(:,:,yr),sim_shocks(:,:,yr),sim_types(:,:,yr),&
            sim_capital(:,:,yr),hbar,w(yr),shocks,types,partprod,&
            r(yr),tr(yr),gov_spend(yr)-Ep(yr)*taue,population(:,2),avg_earnings(yr),ss(:,:,yr),survive,&
            ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))



else if( routine==1  ) then
!rebate
        govt_taxes(yr)=tax_differ2(lambda0(yr),sim_labor(:,:,yr),dummy_zeros,sim_shocks(:,:,yr),sim_types(:,:,yr),sim_capital(:,:,yr),&
            hbar,w(yr),shocks,types,partprod,&
            r(yr),tr(yr),gov_spend(yr),population(:,2),avg_earnings(yr),ss(:,:,yr),lambdak(yr),survive,rebate_level(yr),&
            ss_cap(yr),lambda1(yr))
deallocate(dummy_zeros)
elseif(routine==2) then
    govt_taxes(yr)=0.0D0

else
    print*,"Not Specified Right"
end if

end do
    tax_diff_one=MAXVAL(tax_diff)

        !Check for bad convergence within tax loop

        tax_loop=tax_loop+1
        if (tax_loop > max_tax_loop) then
          tax_diff_one = 0
          bad_combination_inner_nostop=1
        end if
        CALL MPI_Bcast( tax_diff_one, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    END DO
print*,"out of tax loop"

diff_parameters=0.0D0
K1_prime(1)=K_start
tr1(1)=tr_start
K(1)=K_start
tr(1)=tr_start
do yr=1,yr_T
print*,"aggregating",yr
year_calculating=yr
    K1_prime(yr+1)=0.0D0
    N1_prime(yr)=0.0D0
    tr1(yr+1)=0.0D0
    hours_jnr(yr)=0.0D0
    Ec1(yr)=0.0D0

    avg_earnings1(yr)=0.0D0
    avg_earnings_ss1(:,:,yr)=0.0D0


    !Calculate new aggregates
    do it=1,simulationsN
        do s=1,Jnr-1
            avg_earnings1(yr)=avg_earnings1(yr)+sim_earnings(s,it,yr)*population(s,2)/(real(simulationsN)*sum(population(1:Jnr-1,2)))
            if (sim_shocks(Jnr,it,yr) .le. shocksN/2) then
                avg_earnings_ss1(sim_types(1,it,yr),sim_shocks(Jnr,it,yr),yr)=&
                    avg_earnings_ss1(sim_types(1,it,yr),sim_shocks(Jnr,it,yr),yr)+&
                    min(ss_cap(yr),sim_earnings(s,it,yr))*population(s,2)/(ss_count(sim_types(1,it,yr),sim_shocks(Jnr,it,yr),yr)*&
                    sum(population(1:Jnr-1,2)))
            else
                avg_earnings_ss1(sim_types(1,it,yr),sim_shocks(Jnr,it,yr)-shocksN/2,yr)=&
                    avg_earnings_ss1(sim_types(1,it,yr),sim_shocks(Jnr,it,yr)-shocksN/2,yr)+&
                    min(ss_cap(yr),sim_earnings(s,it,yr))*population(s,2)/(ss_count(sim_types(1,it,yr),sim_shocks(Jnr,it,yr)-shocksN/2,yr)*&
                    sum(population(1:Jnr-1,2)))
            end if
            Ec1(yr)=Ec1(yr)+sim_energy(s,it,yr)*population(s,2)/(simulationsN)
        end do
        do s=1,Jnr-1
            if (sim_labor(s,it,yr)<hbar) then
                N1_prime(yr)=N1_prime(yr)+partprod*sim_labor(s,it,yr)*shocks(sim_shocks(s,it,yr))*types(sim_types(s,it,yr))*epsilon1(s,1)*population(s,2)/simulationsN
            else
                N1_prime(yr)=N1_prime(yr)+sim_labor(s,it,yr)*shocks(sim_shocks(s,it,yr))*types(sim_types(s,it,yr))*epsilon1(s,1)*population(s,2)/simulationsN
            end if
            K1_prime(yr+1)=K1_prime(yr+1)+sim_capital(s+1,it,yr+1)*population(s,2)/simulationsN
            tr1(yr+1)=tr1(yr+1)+sim_capital(s+1,it,yr+1)*population(s,2)/simulationsN*(1.0D0-survive(s))
            hours_jnr(yr)=hours_jnr(yr)+sim_labor(s,it,yr)*population(s,2)/simulationsN
        end do
        do s=Jnr,J
            Ec1(yr)=Ec1(yr)+sim_energy(s,it,yr)*population(s,2)/(simulationsN)
            K1_prime(yr+1)=K1_prime(yr+1)+sim_capital(s+1,it,yr+1)*population(s,2)/simulationsN
            tr1(yr+1)=tr1(yr+1)+sim_capital(s+1,it,yr+1)*population(s,2)/simulationsN*(1.0D0-survive(s))
        end do
    end do

print*,"aggregatinga",yr
!Lump sum rebate
if(routine==1) then
    lump_sum(yr)=ind_eng_tax(yr)+taue*E(yr)
elseif (routine==2 .or. routine ==3 .or. routine ==4) then
    lump_sum(yr)=0.0D0
else
    print*,"Not Specified Right"
end if

        do e_it=1,shocksN/2
            do it=1,typesN
                ss1(it,e_it,yr)=ss_start(it,e_it)*(Ce_for_share_start*(taue+pe(yr))+C_for_share_start)/(Ce_for_share_start*pe_start+C_for_share_start)
            end do
        end do




        hours_jnr(yr)=hours_jnr(yr)/sum(population(1:Jnr-1,2))
        K1_prime(yr+1)=K1_prime(yr+1)/(1.0D0+growth)
        tr1(yr+1)=tr1(yr+1)/(1.0D0+growth)

    diff_parameter_v(1,yr)=abs(((tr1(yr+1))-tr(yr+1)))/tr(yr+1)
    diff_parameter_v(2,yr)=abs(((K1_prime(yr+1))-K(yr+1)))/K(yr+1)
    diff_parameter_v(3,yr)=abs(((N1_prime(yr))-N(yr)))/N(yr)
    diff_parameter_v(4,yr)=maxval(abs(((ss1(:,:,yr))-ss(:,:,yr)))/ss(:,:,yr))
    diff_parameter_v(5,yr)=abs(tax_diff_keep(yr))
    diff_parameter_v(6,yr)=max(abs(((avg_earnings1(yr))-avg_earnings(yr)))/avg_earnings(yr),abs(((Ec1(yr))-Ec(yr)))/Ec(yr)&
        ,maxval(abs(((avg_earnings_ss1(:,:,yr))-avg_earnings_ss(:,:,yr)))/avg_earnings_ss1(:,:,yr)))
    diff_parameters(yr)=MAXVAL(diff_parameter_v(:,yr))
    print*,"K",K1_prime(yr),K(yr)
    print*,"K1",K1_prime(yr+1),K(yr+1)
    print*,"N1",N1_prime(yr),N(yr)
    print*,"ec1",Ec1(yr),Ec(yr)
    print*,"tr1",tr1(yr+1),tr(yr+1)
    print*,"ss1",ss1(:,:,yr),ss(:,:,yr)
    print*,"avg_earnings1",avg_earnings1(yr),avg_earnings(yr)
    print*,"hours",hours_jnr(yr)
    print*,"lump sum",lump_sum(yr)
    print*,"diff_parameters_v",diff_parameter_v(:,yr)
print*,"aggregatingb",yr
!!!!Update Aggregates
    Ec(yr)=(Ec1(yr)*a+Ec(yr)*(1.0D0-a))
    N(yr)=(N1_prime(yr)*a+N(yr)*(1.0D0-a))
    K(yr+1)=(K1_prime(yr+1)*a+K(yr+1)*(1.0D0-a))
    ss(:,:,yr)=(ss1(:,:,yr)*a+ss(:,:,yr)*(1.0D0-a))
    tr(yr+1)=(tr1(yr+1)*a+tr(yr+1)*(1.0D0-a))
    avg_earnings(yr)=(avg_earnings1(yr)*a+avg_earnings(yr)*(1-a))
    avg_earnings_ss(:,:,yr)=(avg_earnings_ss1(:,:,yr)*a+avg_earnings_ss(:,:,yr)*(1-a))
print*,"aggregatingc",yr
    if (bad_combination_inner_stop == 1) then
        diff_parameters(yr)=0.0D0
    end if
    taxes_clear(out_loop,yr)=0.0D0
    if (abs(govt_taxes(yr)/gov_spend(yr))>tol) then
        taxes_clear(out_loop,yr)=1.0D0
        non_clear(yr)=1
        diff_parameters(yr)=1.0D0
    end if
    if (taxes_clear(out_loop,yr)==0.0D0) then
        allocate(zeros2(max_outer_loops+1))
        zeros2=0.0D0
        taxes_clear(:,yr)=zeros2
        deallocate(zeros2)
    end if
    if (sum(taxes_clear(:,yr))>max_non_clear) then
        diff_parameters(yr)=0.0D0
        bad_combination_inner_stop=1
    end if
    if (out_loop > max_outer_loops) then
        bad_combination_outer=1
        diff_parameters(yr)=0.0D0
    end if
print*,"aggregatingd",yr
end do
!end if
    diff_parameter=MAXVAL(diff_parameters)
    out_loop=out_loop+1
    print*,"Loop"
    print*,"Plus one"
    print*,rebate_loop,out_loop
    print*,"lambda0",lambda0
    print*,"lambdak",lambdak
    print*,"govt miss",govt_taxes/gov_spend
    print*,"diffparm",diff_parameter
    print*,"taxdiff",tax_diff_keep
    print*,"hours",hours_jnr
    print*,"K",K
    print*,"N",N
    print*,"Y",Y
    print*,"ss",ss(1,:,:)
    print*,"ss",ss(2,:,:)
    print*,"tr",tr
    print*,"E",E
    print*,"Ec",Ec
    print*,"avg_earnings_ss",avg_earnings_ss
    print*,"avg_earnings",avg_earnings
    print*,"w",w
    print*,"r",r

    print*,bad_combination_inner_nostop
    CALL MPI_Bcast( diff_parameter, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )

END DO

do yr=1,yr_T
year_calculating=yr
do it =J,1,-1
    do it2=1,simulationsN
        ypp_temp=sim_earnings(it,it2,yr)-.5D0*tauss(yr)*min(sim_earnings(it,it2,yr),ss_cap(yr))
        sim_labor_tax_new(it,it2,yr)=lab_taxer(ypp_temp,lambda0(yr),lambda1(yr))
        sim_labor_tax_old(it,it2,yr)=old_lab_taxer(ypp_temp)
        if(it<Jnr-1) then
            sim_rebate(it,it2,yr)=rebate_calculator(sim_earnings(it,it2,yr),(sim_capital(it,it2,yr)+tr(yr))*(r(yr)),&
            1,sim_labor_tax_new(it,it2,yr),lambdak(yr)*(sim_capital(it,it2,yr)+tr(yr))*(r(yr)),ss_cap(yr))
        else
            sim_rebate(it,it2,yr)=rebate_calculator(0.0D0,(sim_capital(it,it2,yr)+tr(yr))*(r(yr)),&
            0,0.0D0,lambdak(yr)*(sim_capital(it,it2,yr)+tr(yr))*(r(yr)),ss_cap(yr))
        end if
    end do
end do
end do


!!COmbo



         OPEN(UNIT=12, FILE="/<filepath transition output>/sim_labor_tax_old_combo.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
               do it=1,simulationsN
                do it2=1,yr_T
          WRITE(12,887) sim_labor_tax_old(:,it,it2)
          end do
          end do
          CLOSE(UNIT=12)


         OPEN(UNIT=13, FILE="/<filepath transition output>/sim_labor_tax_new_combo.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
               do it=1,simulationsN
                do it2=1,yr_T
          WRITE(13,887) sim_labor_tax_new(:,it,it2)
          end do
          end do
          CLOSE(UNIT=13)


         OPEN(UNIT=14, FILE="/<filepath transition output>/sim_labor_earnings_combo.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
               do it=1,simulationsN
                do it2=1,yr_T
          WRITE(14,887) sim_earnings(:,it,it2)
          end do
          end do
          CLOSE(UNIT=14)


         OPEN(UNIT=15, FILE="/<filepath transition output>/sim_labor_output_combo.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
               do it=1,simulationsN
                do it2=1,yr_T
          WRITE(15,887) sim_labor(:,it,it2)
          end do
          end do
          CLOSE(UNIT=15)

          OPEN(UNIT=16, FILE="/<filepath transition output>/sim_energy_output_combo.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
               do it=1,simulationsN
                do it2=1,yr_T
          WRITE(16,887) sim_energy(:,it,it2)
          end do
          end do
          CLOSE(UNIT=16)

         OPEN(UNIT=17, FILE="/<filepath transition output>/sim_capital_output_combo.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
               do it=1,simulationsN
                do it2=1,yr_T
          WRITE(17,887) sim_capital(:,it,it2)
          end do
          end do
          CLOSE(UNIT=17)

         OPEN(UNIT=18, FILE="/<filepath transition output>/sim_consumption_output_combo.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
               do it=1,simulationsN
                do it2=1,yr_T
          WRITE(18,887) sim_consumption(:,it,it2)
          end do
          end do
          CLOSE(UNIT=18)


        OPEN(UNIT=21, FILE="/<filepath transition output>/sim_shocks_output_combo.dat",STATUS="replace", &
               FORM="formatted")
               do it=1,simulationsN
                 do it2=1,yr_T
         WRITE(21,886) sim_shocks(:,it,it2)
          end do
          end do
          CLOSE(UNIT=21)

         OPEN(UNIT=18, FILE="/<filepath transition output>/sim_rebate_output_combo.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
               do it=1,simulationsN
                do it2=1,yr_T
          WRITE(18,887) sim_rebate(:,it,it2)
          end do
          end do
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<filepath transition output>/K.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) K
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<filepath transition output>/N.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) N
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<filepath transition output>/Ec.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) Ec
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<filepath transition output>/tr.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) tr
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<filepath transition output>/lambda0.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) lambda0
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<filepath transition output>/lambda1.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) lambda1
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<filepath transition output>/lambdak.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) lambdak
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<filepath transition output>/rebate_slope.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) rebate_slope
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<filepath transition output>/rebate_level.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) rebate_level
          CLOSE(UNIT=18)


         OPEN(UNIT=18, FILE="/<filepath transition output>/avg_earnings.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
          WRITE(18,887) avg_earnings
          CLOSE(UNIT=18)

         OPEN(UNIT=18, FILE="/<filepath transition output>/sim_rebate_output_combo.dat", ACTION="write", STATUS="replace", &
               FORM="formatted")
               do it=1,typesN
                do it2=1,yr_T
          WRITE(18,887) avg_earnings_ss(it,:,it2)
          end do
          end do
          CLOSE(UNIT=18)






do yr=1,yr_T

!!!quintile1
C_end=0.0D0
Ce_end=0.0D0
C_start=0.0D0
Ce_start=0.0D0
do it=1,simulationsN
    do s=1,J
        Ce_end(yr)=Ce_end(yr)+sim_energy(s,it,yr)*population(s,2)/real(simulationsN)
    end do
end do


aggregates(1,yr)=w(yr)
aggregates(2,yr)=r(yr)
aggregates(3,yr)=K(yr)
aggregates(4,yr)=N(yr)
aggregates(5,yr)=tr(yr)
aggregates(6,yr)=0.0D0
aggregates(7,yr)=0.0D0
aggregates(8,yr)=gov_spend(yr)
aggregates(9,yr)=tauss(yr)
aggregates(10,yr)=no_carbon_tax
aggregates(11,yr)=E(yr)
aggregates(12,yr)=lambdak(yr)
aggregates(13,yr)=lambda0(yr)
aggregates(14,yr)=lambda1(yr)
aggregates(15,yr)=avg_earnings(yr)
aggregates(16,yr)=rebate_level(yr)
aggregates(17,yr)=rebate_slope(yr)
aggregates(18,yr)=sum(sim_v(:,yr))/real(simulationsN)
aggregates(19,yr)=0.0D0
aggregates(20,yr)=0.0D0
aggregates(21,yr)=0.0D0
aggregates(22,yr)=0.0D0
aggregates(23,yr)=0.0D0
aggregates(24,yr)=0.0D0
aggregates(25,yr)=0.0D0
aggregates(26,yr)=0.0D0
aggregates(27,yr)=0.0D0
aggregates(28,yr)=0.0D0
aggregates(29,yr)=0.0D0
aggregates(30,yr)=0.0D0
aggregates(31,yr)=0.0D0
aggregates(32,yr)=0.0D0
aggregates(33,yr)=0.0D0
aggregates(34,yr)=0.0D0
aggregates(35,yr)=0.0D0
aggregates(36,yr)=0.0D0
aggregates(37,yr)=0.0D0
aggregates(38,yr)=0.0D0
aggregates(39,yr)=0.0D0
aggregates(40,yr)=0.0D0
aggregates(41,yr)=0.0D0
aggregates(42,yr)=0.0D0
aggregates(43,yr)=0.0D0
aggregates(44,yr)=0.0D0
aggregates(45,yr)=0.0D0
aggregates(46,yr)=0.0D0
aggregates(47,yr)=0.0D0
aggregates(48,yr)=0.0D0
aggregates(49,yr)=0.0D0
aggregates(50,yr)=0.0D0
aggregates(56,yr)=out_loop
aggregates(57,yr)=(Ce_end(yr)+E(yr))*taue
aggregates(58,yr)=hours_jnr(yr)
aggregates(59,yr)=0.0D0
aggregates(60,yr)=0.0D0
aggregates(61,yr)=0.0D0
aggregates(62,yr)=0.0D0
aggregates(63,yr)=0.0D0
aggregates(64,yr)=0.0D0
aggregates(79,yr)=Ec(yr)
aggregates(80,yr)=K1(yr)
aggregates(81,yr)=N1(yr)
aggregates(82,yr)=K2(yr)
aggregates(83,yr)=N2(yr)
aggregates(84,yr)=Ep(yr)
aggregates(85,yr)=Pe(yr)
aggregates(86,yr)=sum(sum(sim_consumption(:,:,yr),2)*population(1:J,2))
aggregates(87,yr)=0.0D0
aggregates(88,yr)=sum(sum(sim_energy(:,:,yr),2)*population(1:J,2))
aggregates(89,yr)=0.0D0
aggregates(90,yr)=0.0D0
aggregates(91,yr)=0.0D0
aggregates(92,yr)=govt_taxes(yr)/gov_spend(yr)
aggregates(93,yr)=maxval(sim_capital(:,:,(yr)))

aggregates(94,yr)=routine
aggregates(95,yr)=labor_capital_rebate_base
aggregates(96,yr)=vary_capital
aggregates(97,yr)=no_increase
aggregates(98,yr)=labor_capital_rebate_base
aggregates(99,yr)=rebate_negative_tax
aggregates(100,yr)=negative_taxes
aggregates(101,yr)=rebate_to_old
aggregates(102,yr)=flat_tax
aggregates(103,yr)=bad_combination_outer
aggregates(104,yr)=0.0D0
aggregates(105,yr)=0.0D0
aggregates(106,yr)=0.0D0
aggregates(107,yr)=0.0D0
aggregates(108,yr)=0.0D0
aggregates(109,yr)=0.0D0
aggregates(110,yr)=0.0D0
aggregates(111,yr)=0.0D0
aggregates(112,yr)=0.0D0
aggregates(113,yr)=0.0D0
aggregates(114,yr)=0.0D0
aggregates(115,yr)=0.0D0
aggregates(116,yr)=0.0D0
aggregates(117,yr)=0.0D0
aggregates(118,yr)=0.0D0
aggregates(123,yr)=0.0D0
aggregates(129,yr)=0.0D0
aggregates(135,yr)=0.0D0

end do

print*,"out",aggregates(12:14,rebate_loop),aggregates(16:17,rebate_loop),aggregates(26:29,rebate_loop)


!end do

!1 lump_sum_rebate
!2 throw_away
!3 cap_tax_rebate
!4 labor_rebate
!if(routine==3) then

         OPEN(UNIT=19, FILE="/<filepath transition output>/aggregates.dat",STATUS="replace", &
               FORM="formatted")
      do it2=1,yr_T
          WRITE(19,987) aggregates(:,it2)
          end do


!!
987 format(135F24.16)


886 format(i9)
887 format(F24.16)
889 format(2F24.16)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Slaves
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!Slaves that are going to work
else if (my_id < (num_procs_run)) then
!do rebate_loop=1,slopeN*levelN*progressivityN+(1+progressivityN)*extra
!!!!!!!!!!!!!!!!!!!!!!!! Initialize loop stuff
tax_diff_one=1
diff_parameter=1.0D0


DO WHILE (diff_parameter>tol)
tax_diff_one=1.0D0
    DO WHILE (tax_diff_one >tol)
    do yr=yr_T,1,-1
    year_calculating=yr
            CALL MPI_Bcast( prices, 11, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            CALL MPI_Bcast( ss(:,:,yr), typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            tr(yr)=prices(1)
            lambda0(yr)=prices(2)
            w(yr)=prices(3)
            r(yr)=prices(4)
            tauss(yr)=prices(5)
            rebate_level(yr)=prices(6)
            ss_cap(yr)=prices(7)
            lump_sum(yr)=prices(8)
            avg_earnings(yr)=prices(9)
            lambdak(yr)=prices(10)
            pe(yr)=prices(11)



!Solve Working
        do s=Jnr-1,2,-1
            CALL MPI_Bcast( v_in(:,:,:), typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            v(:,:,:,s+1,yr+1)=v_in

            call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            call MPI_RECV( end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
            allocate(v_alloc(typesN,shocksN,end_row-start_row+1))
            labor_dummy=0.0D0
            kprime_dummy=0
            v_alloc=0.0D0

            do kit=start_row,end_row
                if(kgrid(kit)<0.0D0) then
                    rate_use=r(yr)/survive(s-1)
                else
                    rate_use=r(yr)
                end if
                do e_it=1,shocksN
                    do i=1,typesN
                        c_upper=c_upper_fnc(epsilon1(s,1),kgrid(kit),kgrid,kgridN,tr(yr),&
                            w(yr)*shocks(e_it)*types(i),rate_use,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))
                        temp=int(MINLOC(c_upper,1,mask=c_upper.gt.0.015D0+ebar))
                        if (c_upper(kgridN)>0.015D0+ebar) then
                            temp = kgridN
                        end if
                        if (c_upper(2)<=0.015D0+ebar) then
                            kprime_dummy=0
                            v_alloc(i,e_it,kit-start_row+1)=-1.0D0/0.0D0
                        else
                            CALL working_valuer(s,kgrid(kit),temp,tr(yr),rate_use,&
                            survive(s),&
                            w(yr)*shocks(e_it)*types(i)*epsilon1(s,1),shock_pi(e_it,:,s,i),&
                            v_in(i,:,1:temp),&
                            v_alloc(i,e_it,kit-start_row+1),kprime_dummy,&
                            labor_dummy,&
                            partprod,hbar,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))
                        end if
                    end do
                end do
            end do
            call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
            call MPI_SEND( v_alloc,typesN*shocksN*(end_row-start_row+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
            deallocate(v_alloc)
        end do
        s=1
        CALL MPI_Bcast( v_in(:,:,:), typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
        v(:,:,:,s+1,yr+1)=v_in

end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Simulating
!    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:,:),typesN*shocksN*kgridN*(J+1-jretired)*(yr_T+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
    call MPI_RECV(end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

    allocate(sim_labor_alloc(J,end_row-start_row+1,yr_T))
    allocate(sim_earnings_alloc(J,end_row-start_row+1,yr_T))
    allocate(sim_capital_alloc(J+1,end_row-start_row+1,yr_T+1))
    allocate(sim_consumption_alloc(J,end_row-start_row+1,yr_T))
    allocate(sim_energy_alloc(J,end_row-start_row+1,yr_T))
    allocate(sim_v_alloc(end_row-start_row+1,yr_T))
    allocate(sim_total_income_alloc(J,end_row-start_row+1,yr_T))
    allocate(sim_after_tax_income_alloc(J,end_row-start_row+1,yr_T))
    !!Initialize

    sim_capital_alloc(1,:,:)=0.0D0
    sim_capital_alloc(:,:,1)=sim_capital_start(:,start_row:end_row)
    sim_labor_alloc=0.0D0
    sim_earnings_alloc=0.0D0
    sim_total_income_alloc=0.0D0
    sim_after_tax_income_alloc=0.0D0


    do yr=1,yr_T
    year_calculating=yr
    do it=start_row,end_row
        !Working
        do s=1,Jnr-1
            if(sim_capital_alloc(s,it-start_row+1,yr)<0.0D0) then
                if(s==1) then
                    rate_use=r(yr)/survive(1)
                else
                    rate_use=r(yr)/survive(s-1)
                end if
            else
                rate_use=r(yr)
            end if
            CALL working_valuer(s,sim_capital_alloc(s,it-start_row+1,yr),kgridN,tr(yr),rate_use,survive(s),&
                    w(yr)*shocks(sim_shocks(s,it,yr))*types(sim_types(s,it,yr))*epsilon1(s,1),&
                    shock_pi(sim_shocks(s,it,yr),:,s,sim_types(s,it,yr)),&
                    v(sim_types(s,it,yr),:,:,s+1,yr+1),dummy,sim_capital_alloc(s+1,it-start_row+1,yr+1),&
                    sim_labor_alloc(s,it-start_row+1,yr),partprod,hbar,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))
            if (s==1) then
                sim_v_alloc(it-start_row+1,yr)=dummy
            end if

            if (sim_labor_alloc(s,it-start_row+1,yr)<hbar) then
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1,yr),&
                    sim_capital_alloc(s,it-start_row+1,yr),sim_capital_alloc(s+1,it-start_row+1,yr+1),tr(yr),rate_use,&
                    w(yr)*shocks(sim_shocks(s,it,yr))*types(sim_types(s,it,yr))*epsilon1(s,1)*partprod,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1,yr),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1,yr)=(con_dummy-sim_consumption_alloc(s,it-start_row+1,yr))/(pe(yr)+taue)
                sim_earnings_alloc(s,it-start_row+1,yr)=w(yr)*shocks(sim_shocks(s,it,yr))*types(sim_types(s,it,yr))*&
                    epsilon1(s,1)*partprod*sim_labor_alloc(s,it-start_row+1,yr)
            else
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1,yr),&
                    sim_capital_alloc(s,it-start_row+1,yr),sim_capital_alloc(s+1,it-start_row+1,yr+1),tr(yr),rate_use,&
                    w(yr)*shocks(sim_shocks(s,it,yr))*types(sim_types(s,it,yr))*epsilon1(s,1),ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1,yr),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1,yr)=(con_dummy-sim_consumption_alloc(s,it-start_row+1,yr))/(pe(yr)+taue)
                sim_earnings_alloc(s,it-start_row+1,yr)=w(yr)*shocks(sim_shocks(s,it,yr))*&
                    types(sim_types(s,it,yr))*epsilon1(s,1)*sim_labor_alloc(s,it-start_row+1,yr)
            end if
            sim_total_income_alloc(s,it-start_row+1,yr)=sim_earnings_alloc(s,it-start_row+1,yr)+&
            (sim_capital_alloc(s,it-start_row+1,yr)+tr(yr))*rate_use
            tax_dummy=all_taxer((sim_capital_alloc(s,it-start_row+1,yr)+tr(yr))*rate_use,&
                sim_earnings_alloc(s,it-start_row+1,yr),1,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))+taue*sim_energy_alloc(s,it-start_row+1,yr)
            sim_after_tax_income_alloc(s,it-start_row+1,yr)=sim_total_income_alloc(s,it-start_row+1,yr)-tax_dummy
        end do



        do s=Jnr,J
            sim_labor_alloc(s,it-start_row+1,yr)=0.0D0
            con=-1.0D0
            if(sim_shocks(Jnr,it,yr) .le. shocksN/2) then
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1,yr),&
                    kgrid(capital_positive:),&
                    tr(yr),r(yr),ss(sim_types(s,it,yr),sim_shocks(s,it,yr),yr),avg_earnings(yr),lambdak(yr))
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))

                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1,yr),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it,yr),sim_shocks(s,it,yr),yr),tr(yr),r(yr),survive(s),&
                    v_retired(sim_types(s,it,yr),sim_shocks(s,it,yr),:,s+1-jretired,yr+1),dummy,&
                    sim_capital_alloc(s+1,it-start_row+1,yr+1),avg_earnings(yr),s,lambdak(yr))
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1,yr),&
                    sim_capital_alloc(s+1,it-start_row+1,yr+1),tr(yr),r(yr),ss(sim_types(s,it,yr),sim_shocks(s,it,yr),yr),&
                    avg_earnings(yr),lambdak(yr))
            else
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1,yr),&
                    kgrid(capital_positive:),&
                    tr(yr),r(yr),ss(sim_types(s,it,yr),sim_shocks(s,it,yr)-shocksN/2,yr),avg_earnings(yr),lambdak(yr))
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))

                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1,yr),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it,yr),sim_shocks(s,it,yr)-shocksN/2,yr),tr(yr),r(yr),survive(s),&
                    v_retired(sim_types(s,it,yr),sim_shocks(s,it,yr),:,s+1-jretired,yr+1),dummy,&
                    sim_capital_alloc(s+1,it-start_row+1,yr+1),avg_earnings(yr),s,lambdak(yr))
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1,yr),&
                    sim_capital_alloc(s+1,it-start_row+1,yr+1),tr(yr),r(yr),ss(sim_types(s,it,yr),&
                    sim_shocks(s,it,yr)-shocksN/2,yr),avg_earnings(yr),lambdak(yr))
            end if

            dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1,yr),con_dummy,s)
            sim_energy_alloc(s,it-start_row+1,yr)=(con_dummy-sim_consumption_alloc(s,it-start_row+1,yr))/(pe(yr)+taue)
            sim_total_income_alloc(s,it-start_row+1,yr)=0.0D0+&
            (sim_capital_alloc(s,it-start_row+1,yr)+tr(yr))*r(yr)
            tax_dummy=all_taxer((sim_capital_alloc(s,it-start_row+1,yr)+tr(yr))*r(yr),0.0D0,0,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))+&
                taue*sim_energy_alloc(s,it-start_row+1,yr)
            sim_after_tax_income_alloc(s,it-start_row+1,yr)=sim_total_income_alloc(s,it-start_row+1,yr)
!            -tax_dummy+&
!                rebate_calculator(0.0D0,(sim_capital_alloc(s,it-start_row+1)+tr)*r,0,0.0D0,tax_dummy)
        end do
    end do
    end do

    call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_labor_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_earnings_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_capital_alloc,(J+1)*(end_row-start_row+1)*(yr_T+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_consumption_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_energy_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_v_alloc,(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_total_income_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_after_tax_income_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)


    deallocate(sim_labor_alloc)
    deallocate(sim_capital_alloc)
    deallocate(sim_energy_alloc)
    deallocate(sim_earnings_alloc)
    deallocate(sim_consumption_alloc)
    deallocate(sim_v_alloc)
    deallocate(sim_after_tax_income_alloc)
    deallocate(sim_total_income_alloc)

    CALL MPI_Bcast( tax_diff_one, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end do
call MPI_Bcast( diff_parameter,1,MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
end do

!end do


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Slaves that are not going to work for some
else if (my_id > (num_procs_run-1)) then


!!!!!!!!!!!!!!!!!!!!!!!! Initialize loop stuff
tax_diff_one=1
diff_parameter=1.0D0


DO WHILE (diff_parameter>tol)
tax_diff_one=1.0D0
    DO WHILE (tax_diff_one >tol)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Solving Value Function (These are not working)
    do yr=yr_T,1,-1
    year_calculating=yr
            CALL MPI_Bcast( prices, 11, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            CALL MPI_Bcast( ss(:,:,yr), typesN*(shocksN/2), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            tr(yr)=prices(1)
            lambda0(yr)=prices(2)
            w(yr)=prices(3)
            r(yr)=prices(4)
            tauss(yr)=prices(5)
            rebate_level(yr)=prices(6)
            ss_cap(yr)=prices(7)
            lump_sum(yr)=prices(8)
            avg_earnings(yr)=prices(9)
            lambdak(yr)=prices(10)
            pe(yr)=prices(11)


        do s=Jnr-1,2,-1
            CALL MPI_Bcast( v_in(:,:,:), typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
            v(:,:,:,s+1,yr+1)=v_in
       end do
        s=1
        CALL MPI_Bcast( v_in(:,:,:), typesN*shocksN*kgridN, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
        v(:,:,:,s+1,yr+1)=v_in

    end do

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Simulating
!    CALL MPI_Bcast( v(:,:,:,:),typesN*shocksN*kgridN*(J+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    CALL MPI_Bcast( v_retired(:,:,:,:,:),typesN*shocksN*kgridN*(J+1-jretired)*(yr_T+1), MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    call MPI_RECV(start_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)
    call MPI_RECV(end_row, 1, MPI_INTEGER, root_process, MPI_ANY_TAG, MPI_COMM_WORLD, status, ierr)

    allocate(sim_labor_alloc(J,end_row-start_row+1,yr_T))
    allocate(sim_earnings_alloc(J,end_row-start_row+1,yr_T))
    allocate(sim_capital_alloc(J+1,end_row-start_row+1,yr_T+1))
    allocate(sim_consumption_alloc(J,end_row-start_row+1,yr_T))
    allocate(sim_energy_alloc(J,end_row-start_row+1,yr_T))
    allocate(sim_v_alloc(end_row-start_row+1,yr_T))
    allocate(sim_total_income_alloc(J,end_row-start_row+1,yr_T))
    allocate(sim_after_tax_income_alloc(J,end_row-start_row+1,yr_T))
    !!Initialize

    sim_capital_alloc(1,:,:)=0.0D0
    sim_capital_alloc(:,:,1)=sim_capital_start(:,start_row:end_row)
    sim_labor_alloc=0.0D0
    sim_earnings_alloc=0.0D0
    sim_total_income_alloc=0.0D0
    sim_after_tax_income_alloc=0.0D0


    do yr=1,yr_T
    year_calculating=yr
    do it=start_row,end_row
        !Working
        do s=1,Jnr-1
            if(sim_capital_alloc(s,it-start_row+1,yr)<0.0D0) then
                if(s==1) then
                    rate_use=r(yr)/survive(1)
                else
                    rate_use=r(yr)/survive(s-1)
                end if
            else
                rate_use=r(yr)
            end if
            CALL working_valuer(s,sim_capital_alloc(s,it-start_row+1,yr),kgridN,tr(yr),rate_use,survive(s),&
                    w(yr)*shocks(sim_shocks(s,it,yr))*types(sim_types(s,it,yr))*epsilon1(s,1),&
                    shock_pi(sim_shocks(s,it,yr),:,s,sim_types(s,it,yr)),&
                    v(sim_types(s,it,yr),:,:,s+1,yr+1),dummy,sim_capital_alloc(s+1,it-start_row+1,yr+1),&
                    sim_labor_alloc(s,it-start_row+1,yr),partprod,hbar,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))
            if (s==1) then
                sim_v_alloc(it-start_row+1,yr)=dummy
            end if

            if (sim_labor_alloc(s,it-start_row+1,yr)<hbar) then
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1,yr),&
                    sim_capital_alloc(s,it-start_row+1,yr),sim_capital_alloc(s+1,it-start_row+1,yr+1),tr(yr),rate_use,&
                    w(yr)*shocks(sim_shocks(s,it,yr))*types(sim_types(s,it,yr))*epsilon1(s,1)*partprod,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1,yr),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1,yr)=(con_dummy-sim_consumption_alloc(s,it-start_row+1,yr))/(pe(yr)+taue)
                sim_earnings_alloc(s,it-start_row+1,yr)=w(yr)*shocks(sim_shocks(s,it,yr))*types(sim_types(s,it,yr))*&
                    epsilon1(s,1)*partprod*sim_labor_alloc(s,it-start_row+1,yr)
            else
                con_dummy=workcon(sim_labor_alloc(s,it-start_row+1,yr),&
                    sim_capital_alloc(s,it-start_row+1,yr),sim_capital_alloc(s+1,it-start_row+1,yr+1),tr(yr),rate_use,&
                    w(yr)*shocks(sim_shocks(s,it,yr))*types(sim_types(s,it,yr))*epsilon1(s,1),ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))
                dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1,yr),con_dummy,s)
                sim_energy_alloc(s,it-start_row+1,yr)=(con_dummy-sim_consumption_alloc(s,it-start_row+1,yr))/(pe(yr)+taue)
                sim_earnings_alloc(s,it-start_row+1,yr)=w(yr)*shocks(sim_shocks(s,it,yr))*&
                    types(sim_types(s,it,yr))*epsilon1(s,1)*sim_labor_alloc(s,it-start_row+1,yr)
            end if
            sim_total_income_alloc(s,it-start_row+1,yr)=sim_earnings_alloc(s,it-start_row+1,yr)+&
            (sim_capital_alloc(s,it-start_row+1,yr)+tr(yr))*rate_use
            tax_dummy=all_taxer((sim_capital_alloc(s,it-start_row+1,yr)+tr(yr))*rate_use,&
                sim_earnings_alloc(s,it-start_row+1,yr),1,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))+taue*sim_energy_alloc(s,it-start_row+1,yr)
            sim_after_tax_income_alloc(s,it-start_row+1,yr)=sim_total_income_alloc(s,it-start_row+1,yr)-tax_dummy
        end do



        do s=Jnr,J
            sim_labor_alloc(s,it-start_row+1,yr)=0.0D0
            con=-1.0D0
            if(sim_shocks(Jnr,it,yr) .le. shocksN/2) then
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1,yr),&
                    kgrid(capital_positive:),&
                    tr(yr),r(yr),ss(sim_types(s,it,yr),sim_shocks(s,it,yr),yr),avg_earnings(yr),lambdak(yr))
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))

                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1,yr),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it,yr),sim_shocks(s,it,yr),yr),tr(yr),r(yr),survive(s),&
                    v_retired(sim_types(s,it,yr),sim_shocks(s,it,yr),:,s+1-jretired,yr+1),dummy,&
                    sim_capital_alloc(s+1,it-start_row+1,yr+1),avg_earnings(yr),s,lambdak(yr))
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1,yr),&
                    sim_capital_alloc(s+1,it-start_row+1,yr+1),tr(yr),r(yr),ss(sim_types(s,it,yr),sim_shocks(s,it,yr),yr),&
                    avg_earnings(yr),lambdak(yr))
            else
                con(capital_positive:)=retcon_ar(kgridN-capital_positive+1,sim_capital_alloc(s,it-start_row+1,yr),&
                    kgrid(capital_positive:),&
                    tr(yr),r(yr),ss(sim_types(s,it,yr),sim_shocks(s,it,yr)-shocksN/2,yr),avg_earnings(yr),lambdak(yr))
                con_pos=int(MINLOC(con(:),1,mask=con(:).gt.0.015+ebar))

                CALL retired_valuer(sim_capital_alloc(s,it-start_row+1,yr),0.0D0,kgrid(con_pos),&
                    ss(sim_types(s,it,yr),sim_shocks(s,it,yr)-shocksN/2,yr),tr(yr),r(yr),survive(s),&
                    v_retired(sim_types(s,it,yr),sim_shocks(s,it,yr),:,s+1-jretired,yr+1),dummy,&
                    sim_capital_alloc(s+1,it-start_row+1,yr+1),avg_earnings(yr),s,lambdak(yr))
                con_dummy=retcon(sim_capital_alloc(s,it-start_row+1,yr),&
                    sim_capital_alloc(s+1,it-start_row+1,yr+1),tr(yr),r(yr),ss(sim_types(s,it,yr),&
                    sim_shocks(s,it,yr)-shocksN/2,yr),avg_earnings(yr),lambdak(yr))
            end if

            dummy=energy_sharer(0.00D0,con_dummy*.15D0,con_dummy,energy_share_equation,.0001D0,&
                    sim_consumption_alloc(s,it-start_row+1,yr),con_dummy,s)
            sim_energy_alloc(s,it-start_row+1,yr)=(con_dummy-sim_consumption_alloc(s,it-start_row+1,yr))/(pe(yr)+taue)
            sim_total_income_alloc(s,it-start_row+1,yr)=0.0D0+&
            (sim_capital_alloc(s,it-start_row+1,yr)+tr(yr))*r(yr)
            tax_dummy=all_taxer((sim_capital_alloc(s,it-start_row+1,yr)+tr(yr))*r(yr),0.0D0,0,ss_cap(yr),lambda0(yr),lambda1(yr),lambdak(yr))+&
                taue*sim_energy_alloc(s,it-start_row+1,yr)
            sim_after_tax_income_alloc(s,it-start_row+1,yr)=sim_total_income_alloc(s,it-start_row+1,yr)
!            -tax_dummy+&
!                rebate_calculator(0.0D0,(sim_capital_alloc(s,it-start_row+1)+tr)*r,0,0.0D0,tax_dummy)
        end do
    end do
    end do

    call MPI_SEND( start_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( end_row,1, MPI_INTEGER, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_labor_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_earnings_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_capital_alloc,(J+1)*(end_row-start_row+1)*(yr_T+1), MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_consumption_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_energy_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_v_alloc,(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_total_income_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)
    call MPI_SEND( sim_after_tax_income_alloc,(J)*(end_row-start_row+1)*yr_T, MPI_DOUBLE_PRECISION, root_process, my_id, MPI_COMM_WORLD, ierr)


    deallocate(sim_labor_alloc)
    deallocate(sim_capital_alloc)
    deallocate(sim_energy_alloc)
    deallocate(sim_earnings_alloc)
    deallocate(sim_consumption_alloc)
    deallocate(sim_v_alloc)
    deallocate(sim_after_tax_income_alloc)
    deallocate(sim_total_income_alloc)

    CALL MPI_Bcast( tax_diff_one, 1, MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
    end do
call MPI_Bcast( diff_parameter,1,MPI_DOUBLE_PRECISION, root_process, MPI_COMM_WORLD,ierr )
end do






end if
deallocate(population)

deallocate(v)
deallocate(v_retired)
deallocate(sim_labor)
deallocate(sim_capital)
deallocate(sim_consumption)
deallocate(sim_energy)
deallocate(sim_rebate)






call MPI_FINALIZE(ierr)

end PROGRAM hetero_ss
